uncertainty NMPC-based control scheme for a semi-batch reactor under parameter Computers and Chemical Engineering

(1)

ContentslistsavailableatScienceDirect

Computers and Chemical Engineering

journalhomepage:www.elsevier.com/locate/compchemeng

NMPC-based control scheme for a semi-batch reactor under parameter uncertainty

Alex Kummer

^∗

, Lajos Nagy , Tamás Varga

Institute of Chemical and Process Engineering, University of Pannonia, Veszprém H-8200, Hungary

a rt i c l e i n f o

Article history:

Received 9 April 2020 Revised 4 June 2020 Accepted 27 June 2020 Available online 30 June 2020 Keywords:

Multi-Stage NMPc Worst case scenario Runaway criterion Process safety Parameter identiﬁcation

a b s t r a c t

ExothermicreactionsareoftenperformedinSBRbecausethegeneratedreactionheatcanbemoreeasily keptundercontrolinsuchconstruction.However,anunsuitablecontrolsystemcanleadtothedevelop- mentofthermal runaway,whichmaycauselethaldamage.NMPCwithimplementedthermalrunaway criteriaisapromisingtooltooperateSBRs.However,engineersshouldalwaysconsiderplant-modelmis- matchbecauseuncertainpredictionscancauseundesirablescenarios.Anovelcontrolframeworkispro- posedtooperateSBRsandconsistsofNMPCwiththeimplementedrunawaycriterion,extendedKalman ﬁlterandparameteridentiﬁcationalgorithm.BothMulti-StageNMPCandNMPCwiththeworst-casesce- narioareinvestigatedandtestedintermsofabilitytohandleparameteruncertainty.Theformeris38 timesslowerthanthelatterwithnonoticeableincreaseinreactorperformance.NMPCinitializedbased onthe worst-casescenario withupdating uncertainkineticparameters results inapromising control structureforSBRs.

ThisisanopenaccessarticleundertheCCBY-NC-NDlicense.

(http://creativecommons.org/licenses/by-nc-nd/4.0/)

1. Introduction

Semi-batch reactors (SBRs) are widely applied for highly exothermic chemical reactions so operators canmore easily keep the reactions under control. If one of the reagents is slowly fed to the other component(s), which is already in the reactor, the feedrateofreagentsaffects theheatevolution, soitcan be controlled, anda suitable coolingsystemcan bedesignedto remove the evolvingreactionheat (WesterterpandMolga, 2004). Forex- ample,theoxidationof2-octanolwithnitricacid(vanWoezikand Westerterp, 2000), Williams-Otto process (Rossi et al., 2016) and synthesis of lithium-etinolate are performed in semi-batch reactors.

The designof thereactor structurewithall necessary compo- nentsandthedesignoftheoperatingstrategyshouldbebasedon themathematicalmodeloftheproposedsystem. Thefeedingtra- jectories are generallyquite simpleusing a constant feeding rate over theentireprocess.Based ontheinitial conditionsandprop- ertiesofthereactorsystem, safeoperatingparameters canbede- ﬁnedwithaconstantfeedingratebycalculatingasafetyboundary diagram(Zhangetal.,2019).Aconstantfeedingratecanresultina

∗ Corresponding author.

E-mail address: kummera@fmt.uni-pannon.hu (A. Kummer).

higherbatchtimethanavaryingfeedingrateduringtheoperation.

However,theconstantfeedingrateisnecessarytokeepthesystem inthesafeoperatingregimeifourinformationabouttherealpro- cessescontainsuncertainties,whichwedonotaddressduringthe operation.

Moreover, engineers should always plan forplant-model mis- match, since it is diﬃcult to obtain a model that describes the plantwithsuﬃcientaccuracy. Theplant-model mismatchcan re- sultinan undesirableeventduringtheoperation,anditspreven- tionisnecessary.Manymethodscanbefound intheliterature to handlethisproblem, suchasconsideringuncertainparameters or applyingstate observers.ModelPredictiveControl(MPC)canpro- videarobust controlapproachtohandleuncertaintiesofthesys- tem,wherethefeedingratecanbeoptimized.MPCisanadvanced control systemand can handlesystem boundaries(Rawlingsand Mayne,2009).AnexcellentreviewonthehistoryofindustrialMPC applications can be found in (Qin and Badgwell, 2003). Param- eter uncertainty can be considered by applying the well-known min-max formulation (Kühl et al., 2007), multi-stage methods (Lucia etal., 2013), or tube-based methods (Mayne etal., 2005).

Min-maxMPCtakesintoaccounttheworst-caserealizationofthe parameteruncertainty,although itis conservativeandmayresult inaninfeasibleoptimizationproblem(ScokaertandMayne,1998).

Theconservativenessofmin-maxMPCwasreducedbytakinginto https://doi.org/10.1016/j.compchemeng.2020.106998

(2)

Nomenclature

c concentration[kmol/m³]

cp heat capacity [Component A…G: kJ/kmolK; water:

kJ/kgK]

d reactordiameter[m]

e controlerror

E activationenergy[kJ/kmolK]

F feedrate[m³/s]

F_j coolantfeedrate[m³/s]

h reactorheight[m]

i ithscenario I runawayindication k kthtimeinstance k₀ pre-exponentialfactor

MAT maximumallowedtemperature[°C]

n moles[kmol]

N numberofscenarios OP controlleroutput

p parameters

PST processsafetytime[s]

PSTc criticalprocesssafetytime[s]

r reactionrate[kmol/m³s]

r_c,r_T partial derivativesofthereactionratewithrespect tocandT

Sp componentfractionofcomponent“P”

t time[s]

T temperature[°C]

u controlinput

U heattransferparameter[kW/m²K]

RD relativedeviation qgen generatedheat[kJ/s]

qrem removedheat[kJ/s]

x statevector

x_KC conversionofkeycomponent V reactorliquidvolume[m³] V_j jacketvolume[m³] w weightfactor

^Hr reactionheat[kJ/kmol]

ρ

^density^[kg/m³^]

ρ

ⁱ ^weight^factor^of^ith^scenarioⁱⁿ^PHA

λ

ⁱ ^weight^factor^of^ith^scenarioⁱⁿ^PHA

Subscripts

0 initialconditions control controlhorizon C coolingagent dos dosing

F feed

dos dosed

i ithcomponent

in,i feedconcentrationofithcomponent

J jacket

pred predictionhorizon

R reactor

account the futurefeedback information (Thangavel et al., 2018).

Multi-stageMPCrealizestheuncertaintybyatreeofdiscretesce- narios,whereeachscenariomustsatisfythepredeﬁnedconstraints (Lucia et al., 2013). Puschkle and Misos proposed a robust feasible multi-stage economic nonlinear model-predictive controller (eNMPC)withaheuristicmulti-modelapproach,wheretheworst- casescenariosaregeneratedbasedonsensitivities.Theyneglected thescenariosonthe edgesoftheuncertaintyset withlow sensi- tivity(Puschke andMitsos, 2018). A review ofeMPC is found in

(Ellisetal.,2014). Holtorfetal.presentedMulti-StageNMPCwith on-linegeneratedscenariotreesthatdonotdirectlyscalewiththe numberofuncertainparameters(Holtorfetal.,2019).

In reactoroperations, the first objective should always be the safe operation, especially when SBRs perform highly exothermic reactions, since a poorly designed reactor or poorly chosen operating parameters can lead to the development of reactor runaway. Although the phenomenon of thermal reactor runaway is well known, lethal accidents unfortunately still occurred in the recent past. The root cause of Seveso-disaster was thermal runaway in1976(Fabianoetal.,2017).In 2001,a polymerizationre- actorexplodedbecauseofareactorrunaway(KaoandHu,2002), and in 2007, an explosion occurred at T2 Laboratories, which causedthedeathoffourpeople(Hall,2010).Thermalrunawayre- sultsin a rapidand significant temperatureincrease inthe reactor, which can lead to explosion through the vaporization of re- actants. The temperature increase in the reactor accelerates the heat generation because of the higher reaction rate, which fur- therincreasesthetemperatureandmayresultinreactorrunaway (Schweitzeretal.,2010).Thereareseveralrunawaycriteria,which classifythe reactionoperationstatesasrunawayornon-runaway, so theseequationscan be applied to predictthe developmentof thermal runaway (Kummer and Varga, 2019). Tailored runaway criteria were developed to obtain a more system specific critical equation(Kummeretal.,2019).

Adequatemodelsofbatchreactorscanbe usedforanonlinear model predictive control (NMPC) (Findeisen et al., 2007). More- over, batch andsemi-batch reactors that perform highlyexother- micreactionscanhavehighlynonlinearprocessdynamics,andthe controllermustcopewithit.NMPCcanbeasuitabletool tohan- dle nonlinear processesand is gainingmore attentionbecause it can capture detailed nonlinear dynamics of the system through- outtheentirestatespace(Sekietal.,2001;YuandBiegler,2019).

Since the goal of reactor operation is to maximize productivity while keeping the operation safe in the entire production time, preventionofthermalrunawayisnecessary.Thermalrunawaycri- teria(Strozzi-Zaldivar andModiﬁedDynamicCondition) wereimplemented in NMPCto reliably indicatethe developmentof runaway(Kummeretal.,2020).

Differentstabilityanalysestopredictthedevelopmentofther- mal runaway were successfully implemented in NMPC, such as the batch simultaneous model-based optimization and control (BSMBO&C)algorithm.ThisalgorithmisanextensionofNMPCand dynamic real-time optimization (DRTO) techniques, which use a Booleantermthatpenalizestheobjectivefunctionwhenthecon- trollersystemisclosetothermalrunaway(Rossietal.,2015).Spe- ciﬁcclassesofdeterministic NMPC/DRTOframeworkscanidentify reactor runaways under parameter uncertainty too (Rossi et al., 2017). Strozzi-Zaldivarcriterion can be toostrict; hence, itis not suitable to analyse the stability of semi-batch reactors in some cases (Kummer and Varga, 2017). Kähm-Vassiliadis criterion for exothermicbatch reactors wasintroduced toovercomethisprob- lem, and the proposed stability criterion can be successfully ap- pliedin batchreactor control toperformhighlyexothermic reactions(KähmandVassiliadis,2018a).Thismethodwasalsosuccess- fully generalized to calculate withmultiple reactions (Kähm and Vassiliadis, 2019). Their stability criterion was applied to an in- dustrialcasestudyandtheyconsideredtheparameteruncertainty duringthe process control (Kanavalau et al., 2019). Lyapunovex- ponents asan indicator of stability were successfully realized in NMPC to control batch reactors (Kähm and Vassiliadis, 2018b).

The operation of an industrial semi-batch polymerization reactor wasoptimizedbyconsideringacoolingsystemfailure (Abeletal., 2000). The interaction between control and safety systems was also studied, where an LMPC (Lyapunov-based MPC) systemwas integratedwiththeactivationofasafetysysteminaCSTRtoavoid

(3)

thermal runaway (Zhang et al., 2018), and a good review about LMPC can be found in (Albalawiet al., 2018). Recently, two new NMPC-basedmethodswereintroducedtosolvetheclosed-loopdy- namicoptimizationproblems,whichweretestedonasemi-batch reactorwithpotentialrunawayreactions,wheretheadiabatictem- perature rise was considered to avoid reactor runaway. The ﬁrst methodisbasedonan adaptivebackingoff oftheir boundsalong the moving horizon with a decreasing degree of severity. The second method is a chance-constrained control approach, which considers the relation betweenthe uncertain input and the constrained output variables. Bothmethods considerthe unexpected disturbancesinadvance,whichresultsinarobustcontrolapproach (Arellano-Garciaetal.,2020).

We proposea controlstrategy forSBRs thatperform potential runawayreactions,wheretheparameteruncertaintyisconsidered.

Weconsiderthemultiplicativeuncertainty,sothemodelmatrices areuncertain.Themostcrucialuncertaintysourcesfromthereac- torrunawayviewpointcan bethe kineticparameters, heattrans- fer parametersandmixingeﬃciencyforSBRs. Theuncertaintyin kinetic parameters is considered the source of model-plant mis- matchin ourwork. NMPC withtheembedded runawaycriterion can be an effective tool to keep the reactor in the controllable zone,whiletheproductivityismaximized.Weinvestigatedworst- case scenario andMulti-Stage NMPC to handle parameter uncertainty,althoughthecomputationcostofMS-NMPCismuchhigher.

ThegoalistodevelopacontrolframeworkforSBRswithexother- micreactions,whichcanbeappliedonlineinrealreactorsystems, sothecomputationtimeiscritical.Therefore,weinvestigatedthe worst-casescenariowithiterativelyupdatinguncertainparameters bytheleast-squaresmethod.NMPCnaturallyincludesthemodelof theprocess,althoughtherealprocessinthiscaseisalsoamodel.

Our future work will be aboutimplementing the proposed con- trolschemeintoareallaboratoryreactorsystem.Currently,wecan generateaplant-modelmismatchtoinvestigatetheproposedcon- trolframework.

2. ProposedcontrolstructureofSBRs

This section will introduce the proposed general control methodology for SBRs that perform potential runaway reactions.

Thepracticalapplicationofthiscontrolstructurewillbepresented inSections3–4.

OurperspectiveisthatusinganNMPCwithimplementedrun- awaycriterioncanbeanexcellentsolutiontocontrolSBRs,butwe must handle parameter uncertainty to develop a reliable indica- tiontool.TheproposedcontrolschemeforSBRsisshowninFig.1. When exothermicreactions occur inthereactor, thepoorly iden- tiﬁed kinetic parameters (or other model parameters) can easily leadtothedevelopmentofthermalrunawayintheapplicationof themostappropriatecriteria.Intheproposedcontrolscheme,the parameter uncertaintyishandledby thecombinationofstate estimation andamodelidentiﬁcationalgorithm.Sincethermalrun- awaycanhavelethalconsequence,theparameteruncertaintymust behandledsothattheprobabilityofrunawayisclosetozero.For thispurpose,Multi-StageNMPCorworst-casescenariocanbeap- plied.Although we investigatedMulti-Stage NMPC, itshigh com- putational cost is not encouraging (see Section 4.2.1);hence, we suggest applyingtheworst-casescenariowithupdatinguncertain parameters(seeSections4.2.2–4.2.3).

InFig.1,uisthecontrolinputs, yisthereactormeasurement outputs, x¯is the estimated states ofthe reactor, andx¯_iis there- quired estimated states for model parameter identiﬁcation. This schemeisa generalrepresentationofthe proposedcontrol structure.Inourcase,uconsistsofOPF andOPC (valvepositionsinthe feedlineandcoolingagentline,respectively).yincludesT_R,T_J,and V_L (reactor temperature, jackettemperatureandliquid volumein

Fig. 1. Proposed control scheme for SBRs.

thereactor).x¯consistsofc¯,T¯R, andT¯J(estimatedconcentration,re- actorandjackettemperature,andreactionrateconstants).x¯_iinour caseincludesk¯(reactionrateconstants),andpconsistsofk¯₀andE¯, whichareidentiﬁed parameters.Inthefollowing sections,we in- troducepartsoftheproposedcontrolstructureinmoredetail.

2.1. TemperaturecontrolofSBR

ThemaingoaloftheSBRoperationissimilartoanyotherkind of reactor orprocess unit, such asto maximize the productivity whilekeepingthereactorsafeduringtheentireoperation.There- actor temperatureis controlled by manipulating the flow rate of feedingreagents(TV001) andcooling agent(TV002)intothereac- torjacket,whichistheoftenappliedscheme.Thecontrolstructure ofthereactorsystemisshowninFig.2.Themass-flowofreagent feedis integrated, theamount offed reagentis calculated (V_dos), andtheliquid levelismeasured withLIT001. Whenthe required reagentisadded(V_dos,0),thereactortemperaturecontrol switches the actuator, and the temperature is controlled by manipulating theflow ofcoolingagentinthemixingphase,which was80% of thetotalcoolingcapacityinthefirstphaseofoperation.

Actuator=

TV001 V_dos<V_dos_,₀

TV002 Vdos,0=Vdos (1) A safe operation is characterized by applying two safety reg- ulations. One ofthem is a predeﬁned Maximum AllowableTem- perature(MAT),sothereactortemperature(TIT001)cannotexceed MAT.Theothersafetyboundisapplyingathermalrunawaycrite- riontoavoid dangerousrunaway states.Runawaycriteriaclassify runawayandnon-runawaystatesbasedonthestatevariablesand parametersofthestudiedprocess(concentrationofreagents,pro- cesstemperature, heat ofreaction, heat transfer parameter, etc.).

Therefore,avoidingrunawaystatesincreasesthesafetyofthereac- toroperation.WeappliedtheModiﬁedDynamicCondition(MDC) in the proposed control scheme (Kummer andVarga, 2019). The MDCstatesthat theproductoftheincreaseinheatremovalwith increasing temperature and the generated to removed heat ratio must exceed the difference between the rate at which the heat generationincreasesasaresultofanincreaseintemperatureand the reduction in the reaction rate as a result of a decrease in reagentconcentration.IftheMDCissatisﬁed.thereactorstatesare considerednon-runaway(i.e.,normal)operation.

∂

^q^gen

∂

^T

c

+

∂

^m

∂

^c

T

≤dqrem

dT qgen

qrem

(2)

(4)

Fig. 2. Proposed control structure of SBRs.

whereq_gen is the generated heat in the reaction, q_rem is the removedheat,andmisamassbalancefunction.

LetusintroducevariableI(k)to evaluatethereactoroperation basedonrunawayandnon-runaway states.IfEq.(2). issatisﬁed, then I(^k)=0 (normal operation); otherwise, I(^k)=1 (runaway), wherekisthek-thtimeinstance.

2.2.Open-loopoptimizationproblem

The goal is to maximize the productivity while thermal runaway doesnot develop, so the conversionof the key component (x_KC)andselectivityfortheproduct(S_P) areconsideredintheob- jectivefunctionnexttotherunawaystates(I_k),andhigherreactor temperaturesthanMAT(e⁺)shouldbeavoided.Theobjectivefunc- tion(orstagecostifwerefertoMulti-StageNMPC)isdenotedby L, whichrepresenta generalcost function.The terms ofthecost functionare weighted (wx,ws, wu, w_I,w_T),so a well performing controlcanbereachedindifferentapplications.Inthethirdterm inEq.(4),signiﬁcantchangesinthemanipulatedvariablesarepe- nalized.

e⁺=max

TR,k−MAT;0

(3) L= −wxx_KC−wsSp+wu

|

^u^k⁻^u^k−1

|

⁺^wÎÎ^k⁺^w^Tê⁺ ⁽⁴⁾

min uk

tpred

0

L

(

^x^k^,^u^k

)

⁽⁵⁾

whichissubjectto

xk+1=f

(

^xk,uk,dk

)

⁽⁶⁾

0≤u_k≤100% (7)

where u_k∈Rⁿ^uu_k is the control input (control valve), nu is the numberof control inputs, each state (x_k₊₁∈Rⁿ^x) is a function of thepreviousstate(x_k),andnxisthenumberofstates.Therealiza- tionofuncertaintyisdenotedasd_k∈Rⁿ^d,wheren_d isthedimen- sionoftheuncertaintyvector.

2.3. NMPCtohandleparameteruncertainty

Wealways mustcount on plant-modelmismatch, sowe must addresstheparameteruncertainty.Thissectionintroducesthefor- mulationofMulti-StageNMPCandWorst-caseScenariotohandle thisproblem.

2.3.1. Multi-StageNMPC

Fig. 2 illustrates how the Multi-Stage NMPC works with the given horizon lengths. For Multi-Stage NMPC, combinations of maximal,minimalandnominalvaluesofuncertainparametersare considered,whichusuallyresultsinarobustbehaviourofthecon- troller(Thangaveletal.,2018).

Eachpathofthescenariotreeiscalledascenarioandindicated asi,anditcontainsallstatesx_k^j andcontrol inputsu_k^j ofscenario i.Thesetofalloccurringindices(j,k)isdenotedbyTLucia,2014).

ThenumberofscenariosisintroducedbyN.Thecostofeachsce- narioisconsideredwiththesameweight,sothemeanvalueofthe costswill givetheobjectivevalue.TheformulationofMulti-Stage NMPCisshowninEqs.(8)–((11).

min u_k^j

N

i=1

1 N

tpred

0

Li

xⁱ_k,uⁱ_k,d_kⁱ

(8)

whichissubjectto x_k^j₊₁= f

x^p_k⁽^j⁾,u_k^j,d_k^j

∀ (

^j,^k⁺¹

)

^∈^T ⁽⁹⁾

0≤u_k^j≤100%

∀ (

^j,^k

)

^∈^T ⁽¹⁰⁾

u_k^j=u^l_ki f x_k^p⁽^j⁾=x_k^p⁽^l⁾

∀ (

^j,^k

)

^,

(

^l,^k

)

^∈^T ⁽¹¹⁾

wherex_k^p⁽ⁱ⁾istheparentnode.Tocorrectlyrepresentthereal-time decisionproblem,thecontrolinputscannotanticipatethevaluesof theuncertaintythat are realizedafterthecorresponding decision point.Itisimportantbecauseitisnotpossibletogivemultiplein- putvariationstotheprocessatthecurrentstate(PuschkeandMit- sos,2018).ThisconditionisenforcedbyEq.(11),whichrepresents

(5)

Fig. 3. Tree representation of the uncertainty evolution for a Multi-Stage NMPC ( Thangavel et al., 2018 ).

thenon-anticipativityconstraintsthatrequireallcontrolinputsat the same node to be equal. In Fig.3, thiscondition impliesthat u¹₀=u²₀=u³₀; u¹₁=u²₁=u³₁;... [36].

Theoptimizationproblemwassolvedbythemodifiedprogres- sivehedgingalgorithm,whichisadecompositionalgorithm,where non-anticipativity constraints are relaxed by penalizing the difference between the control inputs that should satisfy the non- anticipativity constraints. Its advantage is that the scenarios can be independently solved, so the following Eqs. (12)–((15)) optimization problem must be solved. As shown by S. Lucia applying alonger robusthorizonofone doesnot significantly improve the result,butit requiresmorecomputational effort, becausethe size of the optimization problem increases exponentially with it (Lucia,2014).Sincethelengthofrobusthorizonisoneinthiscase, onlythefirstcontrolinputs(uⁱ₀)ofdifferentscenariosmustsatisfy thenon-anticipativityconstraint.

min u_k^jL_i

xⁱ_k,uⁱ_k,dⁱ_k

+

λ

ⁱ

uⁱ₀−uˆⁱ₀

+

ρ

ⁱ

uⁱ₀−uˆⁱ₀

2

(12)

whichissubjectto x_k^j₊₁= f

x_k^p⁽^j⁾,u_k^j,d_k^j

∀ (

^j,^k⁺¹

)

^∈^T ⁽¹³⁾

0≤u_k^j≤100%

∀ (

^j,^k

)

^∈^T ⁽¹⁴⁾

ˆ uⁱ₀=

N

i=1

1

Nuⁱ₀ (15)

where uˆⁱ₀ is the ﬁctious value towards which the control inputs converge to satisfy the anticipativity constraints. Parameters

λ

ⁱ

and

ρ

ⁱ âreûpdated âtêachîteration ^toîmprove^theconvergence, wheretheupdateruleis:

λ

ⁱ⁼

λ

ⁱ⁺

ρ

ⁱ

uⁱ₀−uˆⁱ₀

(16)

ρ

ⁱ⁼^min

βρ

ⁱ^,

ρ

^max

(17) where

β

^determines ^the ^increase ^of

ρ

ⁱ^. Êqs. ⁽¹²⁾^–(15) âre îter-

ativelysolved until max(^uⁱ1−uˆⁱ₁)<

ε

^. ^After ^several^iterations ^the

non-anticipativity constraints are satisﬁed with desiredtolerance

ε

^.

2.3.2. Worst-casescenario

Two non-desired scenarios can be distinguished. In the first scenario,thereactionrateismuchhigherthanexpected;then,the generatedheatwillbemuchhigherintheprocessthanthemodel, whichmaycausethermalrunaway.Inthesecondscenario,thefed reagentaccumulates because the reaction rateis lower than expected. When the concentration increases to a critical point, the reaction ignites, and thermal runaway may occur. To select the worst case scenario. we must choose between these two possi- blescenarios.Thereisahugedifference inignitiontimebetween thesetwo scenarios. The first scenario results in a more conser- vativesolution,sincethe ignitiontime islower; hence,ifwe can handlethe first scenario,we can alsoavoid the second scenario.

Therefore,we suggest kinetic parametersfor worstcases that re- sultinhigherreaction rates,andwe will apply theﬁrst scenario astheworstcasebyincreasingthepre-exponentialfactor(k₀)and decreasingtheactivation energy(E)totheedgeoftheconﬁdence interval.

2.4.Parameteridentiﬁcation

The logarithm ofthe reaction rate constant (k) linearlyvaries withthe reciprocaloftemperature(Eq.(18)), sotheleastsquares methodcanbeappliedtoestimate thepre-exponential factorand activationenergyofthereaction.

ln

_¯

k

=ln

_¯

k0

− E¯

RT (18)

Estimated reaction rate constants are required to calculate Eq.(18),sothestateobserverisnecessaryinthecontrolstructure.

Since the reliability ofthe estimated kinetic parameters depends onthereliabilityofthestate observer,weimplementacondition that mustbe satisﬁedtooverwritethe actual kineticparameters.

RD=100

σ

^p

μ

^p ^<^1% ⁽¹⁹⁾

where

μ

^p ^and

σ

^p âre ^the^meanând^standard^deviationôfûncer-

tain parameters. Due to thiscondition, the estimatedkinetic parameters do not signiﬁcantly vary, so our reliability in thesepa- rametersishigher.Iftheestimatedparametersare farawayfrom theworstcasescenariotheuncertainparametersarenotupdated.

Therefore,weupdate theseparameters ifEq.(19)issatisﬁedand iftheestimatedvaluesarewithintheconﬁdenceinterval.

2.5.StateestimationbasedontheextendedKalmanﬁlter

In realsystems only some measurements are available online (see Section 2.1), andusually eachor noneofthe concentrations cannotbemeasuredonline.Therefore,stateestimationofthesys- tem is necessary to use an effective NMPC in real systems. We mustestimatetheconcentrationofreagentsandproductsanduse it as a feedback in NMPC. The state estimation is necessary to

(6)

identifyuncertainkinetic parameters. The extendedKalman ﬁlter is a suitable algorithm to estimate states of non-linear systems (QuandHahn,2009;KhodadadiandJazayeri-Rad,2011),soweim- plementedthisalgorithm.Ifthereisa closed-formexpressionfor thepredictedstateasafunctionofthepreviousstate(^x^ˆk)^,^controls (u_k)andnoise(w_k),thepredictedstateiscalculatedbyEq.(20). ˆ

x_k₊₁=f

ˆ

x_k,u_k,w_k

(20) Themeasurementisthefunctionofthestate(x_k)andmeasure- mentnoise(v_k).

zk=h

(

^x^k^,

v

k

)

⁽²¹⁾

TypicallyforSBRs,themeasurementvector(z_k)consistsoftem- peratureandlevelmeasurements, andthe stateestimationvector ˆ

x_kconsists of concentrations and temperatures.

To improve the state estimation accuracy, additional variables wereimplemented intothemodelofEKF,which arethereaction rateconstantsofthereactions(k_l).Estimatingnewstatevariables inEKFisnottimeconsuming,anditdoesnotincreasethecompu- tationaltime.Tosolvethisissue,wemustknowhowthereaction rateparametersvarywithtime;sincethereactionrateparameter varieswithtemperature(k = k(T)),the following equation isde- ﬁned:

dk_l dt =dk_l

dT dT

dt =k₀_,lexp

−E_l RT

_E

l

RT² dT

dt (22)

2.6.Processsafetytimetodeﬁnethelengthofthepredictionhorizon

IncaseofMPCthepredictionhorizonmustbesufficientlylong tocapturethefirstrunawaystates,sothesestatescanbeavoided duringthe reactor operation. Runaway criteriacan be applied to definethe length ofpredictionhorizonby defining the firsttime instantwhenarunawaystateoccurs(alsocalledtheprocesssafety time,PST)accordingtotheappliedrunawaycriterion.ThePSTcan bedefinedasfollows:

PST=argmink

(

^I^k⁼¹

)

⁽²³⁾

Since the PSTis a function of state variables and systempa- rameters,we mustcalculatethe PSTfortheworst case, which is calledthecriticalPST.Thus,theSBRsystemmustbeconsidereda batch(i.e.,ifundesiredeventsoccurinSBRs,ourﬁrstactionshould be closing the feeding valve of reagent), and we must calculate withthose state variables where the probability ofaccumulation isthehighest(lowreactortemperature,initialreagentsconcentra- tion,etc.).Differentscenarios(i=1…n)mustbeanalysedbyvary- ingtheinitialconcentrationofreagents(whicharefedintothere- actor)todeﬁnemaximumprocesstemperatures(T⁽ⁱ⁾max)andpro- cesssafetytimes(PST⁽ⁱ⁾). Criticalinitial states(x_0,c) willbe these initialstateswherethemaximumprocesstemperatureequalsMAT.

LengthofpredictionhorizonwillbethePSTatcriticalinitialstates (Kummeretal.,2020).

t_pred= PST

|

i(TR,max=MAT) (24) 3. Processmodelandanalysis

Thissectionpresentstheprocessmodeloftheinvestigatedfed- batch reactorfromWilliams-Otto process, where normalandab- normaloperationsthatcausethermalrunawayarepresented.The processsafetytime ofthe systemisalso calculatedto deﬁne the lengthofthepredictionhorizon.

3.1.Processmodelofthereactorsystem

The Williams-Otto process (WOP) has been used for years to test different control and optimization algorithms Arellano- Garcia et al., 2020). We optimize the fed-batch version of this

process as presentedin (Rossi etal., 2016). In the Williams-Otto process three exothermicreactions occur, which are presentedin Eqs.(25)–((27)followedbytheequationofreactionrates.

A+B→C r1=k0,1exp

−E1

RT

cAcB (25)

B+C→P+E r2=k0,2exp

−E2

RT

cBcC (26)

C+P→G r3=k0,3exp

−E3

RT

cCcP (27)

ComponentAispreloadedandcomponentBiscontinuouslyfed intothereactor.ThedesiredproductiscomponentP,andtwoco- productscanbeformed:componentsEandG.

The following differential equations Eqs. (28)–((31)) describe thedynamicalbehaviourofthereactorsystem:

dc_i dt = Fⁱⁿ

VR

cⁱⁿ_i −ci

+

NR

l=1

υ

ilR_l i=1...NC (28)

dVR

dt =Fⁱⁿ (29)

dTR

dt = 4U DR NC

i=1cicpi

T_j−TR

+Fⁱⁿ ^N_i₌₁^C c^IN_i cpi

VR NC i=1cicpi

Tⁱⁿ−TR

−

NR l=1

^Hr,lR_l

NC i=1cicpi

(30)

dTj

dt = 4UVR

D_RV_j

ρ

jcp_j

T_R−T_j

+ Fj

V_j

T_jⁱⁿ−T_j

(31)

The kinetic parameters, component properties and reactor constructional and operating parameters are summarized in Tables1–3.Theparameterswillbehandledasnominalhereinafter.

TheconstraintsaredeﬁnedinTable3,suchastheMAT,andmax- imumfeedratesofthereagentandcoolingagent.

Table 1

Kinetic and thermodynamic parameters of reactions ( Rossi et al., 2016 ; Sriram and Stevens, 1973 ).

Parameter Value Unit

Pre-exponential factors k 0,1 1.3833 10 ⁵ _kmol^m³_s k 0,2 6.0098 10 ⁷

k 0,3 2.2288 10 ¹¹

Activation energies ^E_R¹ 6450 K

E2

R 8778.5

E3

R 11,155

Heat of reactions H r,1 −1.851 10 ⁵ _kmol^kJ H r,2 −2.5765 10 ⁵

H r,3 −5.053 10 ⁵

Table 2

Component properties ( Rossi et al., 2016 ; Sriram and Stevens, 1973 ).

Component Molecular weight [kg/kmol]

Speciﬁc heat [kJ/kmolK]

A 142 321.204

B 60 127.14

C 202 352.288

E 81 166.212

P 181 426.617

G 383 844.132

[kJ/kgK]

Cooling agent (water) 18 4.186

(7)

Table 3

Reactor constructional and operating parameters ( Rossi et al., 2016 ; Sriram and Stevens, 1973 ).

Parameter Value Unit

d 1 M

h 3.5 M

V j 0.8236 m ³

U 0.8 ^kW_m2

c in 1 ^kmol_m3

T in,R 298 K

T in,j 298 K

c 0,A 1 ^kmol_m3

V 0 0.5 m ³

T 0,R 312 K

T 0,J 308 K

F max 1e-3 ^m_s³

F j,max 1e-2 ^m_s³

MAT 335 K

Fig. 4. Reactor behaviour at different dosing times and cooling agent ﬂow rates.

3.2. AnalysisofwopinSBR

Thebehaviour oftheinvestigatedreactorsystemwasanalysed withnocontrolsystem(i.e.,theBreagentfeedisconstant),where the maximum process temperatures are analysed in functionsof thedosingtimeandﬂowrateofcoolingagent.Theremainingap- plied parameters are shown in Section 3.1. As shown in Fig. 4, poorlychosenoperatingparameterscandevelopthermalrunaway.

Accordingtothemodel,theprocesstemperaturecanexceed900K.

Themaximumtemperaturerapidlyincreases,andthereisnointe- riorpointbetweennormalprocesstemperatures(underMAT)and runawaytemperatures(>900K).Althoughtheoptimalfeedingtra- jectory can increase the productivity, increasing the ﬂow rateof coolingagentenablesanoperationwithlessdosingtime,asshown inFig.4.

RunawaystatesaredistinguishedbytheMDCcriterion,andthe derived criticalequation forthe processisintroduced inEq.(32). nr

l=1

−

^Hr,jr_j,TV− nr

l=1

nc

i=1

r_l_,_c_li≤

nrl=1−

^Hr,lr_lV

TR−T_j (32)

wherenristhenumberofreactions,ncisthenumberofreagents inthe L -threaction,r_T andrc are thederivativesofthereaction

Fig. 5. Reaction rates during an operation (Dosing time: 1.8 hr, Feed rate: 0.55 m ³/hr, Cooling ﬂow rate: 36 m ³/hr).

ratewithrespecttotemperatureandconcentrationofreagentsre- spectively.

To avoid thermal runaway uncertain kinetic parameters are quitesigniﬁcant,soweinvestigatehowtheparallelreactionsdom- inateduringthereactoroperation.Fig.5showsthereactionrates;

thefirstreaction(R1)hasthehighestrateduringthewholeoper- ation.Toinvestigatetheproposedcontrolscheme,weonlychoose thekinetic parameters ofthe first reaction asuncertain. Because thefirstreactionisdominant,theuncertaintyofthisreactionhas thehighesteffectonthebehaviourofthereactor.

3.3.Processsafetytimeofthesystem

Aspresented inSection 2.6,the length ofthepredictionhori- zoncan be deﬁned basedon the process safetytime ofthe system.MaximumreactortemperaturesandPSTsareinvestigated,as showninFig.6.MATisreachedat~1.36 ^kmol_m₃ initialconcentration, wherethePSTis0.59h. Inthiscase,theminimumlength ofthe predictionhorizonis0.59h.

3.4.Stateestimationoftheinvestigatedsystem

Thissection presentstheefficiency ofEKFonthe investigated model system. First Eqs. (28)–(31) were applied to estimate the statesof the system, and the measured variables are the reactor temperature, jackettemperatureand reactionvolume (the inflow rateismeasured, which isthe only parameterthat increases the reaction volume). The results are generated next to 5% parame- terdeviationinpre-exponentialfactorandactivationenergyofthe firstreaction.

As shownin Fig. 7a,the estimations of reagentconcentration are quite poor, whichcan resultin false runaway indicationand thermalrunawayofthesystem.Iftheﬁrstreactionrateparameter (k₁) is estimated the accuracy can be increased, as presented in Fig.7b.Thestateestimationsareacceptablewiththismodiﬁcation.

4. ResultsofNMPC

This section provides the results of NMPC with and without parameteruncertainty. Theperformance improvement duetopa- rameteridentiﬁcationispresented.Theoptimizationproblemwas

(8)

Fig. 6. PST of the system according to the MDC criterion.

solved by the interior-point algorithm, where the algorithm pro- ceedsamovinghorizon(VanderbeiandShanno,1999).Theapplied parameters,whichwere heuristicallyselected, aresummarized in Table4.

Table 4

Parameters of NMPC.

Sample time T 0 100 s

Prediction horizon t pred 2200 s

Control horizon t contr 500 s

Weight factor in Eq. (4) w e 500 Weight factor in Eq. (4) w u 0.01 Weight factor in Eq. (4) w I 100

Fig. 8. Reactor operation with nominal NMPC.

4.1. Resultsoftheopen-loopcontrolwithoutparameteruncertainty

NMPC is tested without anyuncertain parameter and the re- sultsareshowninFig.8.Thereactortemperaturestaysfarbelow

Fig. 7. State estimation based on EKF.

(9)

Fig. 9. Result of MS-NMPC with nominal kinetic parameters.

Fig. 10. Results of NMPC with respect to the worst case.

MATsincetheappliedMDCcriterionconstraintsthereactoroper- ation thatincreasestheprocess safety.At2.5htheconversionof componentAis76%,andtheyieldofPis37%.The averagecom- putationaltimeis11.5speriterationsinthiscase.

4.2. Resultsoftheopen-loopcontrolunderparameteruncertainty

Theeffectoftheparameteruncertaintywasanalysedusingtwo different algorithms inSections 4.2.1 and4.2.2. In the first case, Multi-Stage NMPC was applied, in second case, the worst-case scenario was used to solve the optimization problem under parameteruncertainty. Kineticparameters ofthefirst reaction(k_0,₁, E₁) were chosen as uncertain, where the confidence interval is

±5%.Section 4.2.3 provides theresults of theoptimization prob- lemwhentheuncertainparametersareupdatediteratively.

4.2.1. ResultsofMulti-StageNMPC

The Multi-StageNMPC algorithmwastestedwherethe uncer- tainkineticparameterswerechangedby5%.Twouncertainparam- etersleadtoninescenarios.AsshowninFig.9,thefeedratesare maintainedatlow valuesduetotheuncertainkineticparameters.

The reason is that the constraintsmust be satisﬁed in each scenario,sothe developmentofthermalrunawayisavoided ineach scenario.Therefore,theresultswithMulti-StageNMPCareconser- vativecomparedtothenominalsolution(Fig.8).At2.5hthecon-

Fig. 11. Results of NMPC initialized from the worst-case scenario with updating kinetic parameters.

Fig. 12. Result of the parameter ﬁtness.

version ofcomponent Ais 15.3%, andthe yieldof Pis 2.4%. The averagecomputationtimeis660speriterations,soreal-timeop- timizationisnotfeasiblewiththeMulti-StageNMPCalgorithm.

4.2.2. Resultsoftheworst-casescenario

The worst case is that the real reaction rate is higher than expected, so in the worst-case scenario, the uncertain pre- exponentialfactor increasesby 5% (k_0,1 + 5%),andthe uncertain activationenergy decreases by 5%(E1–5%).The NMPCresultsare showninFig.10,whichnaturallyisaconservativeresult.Thecon- version at2.5h is45%, andtheyield ofproductP is 16%.In the worst-casescenario,theaveragecomputationtimeis17.2sperit- erations,soreal-timeoptimizationisfeasiblewiththisalgorithm.

4.2.3. Worst-casescenariowithupdatinguncertainparameters The results of the proposed control structure, which wasini- tializedfrom theworst case, asshownin Fig. 11. Withupdating uncertainkinetic parameters, thereactortemperature controlbe- comeslessconservativeandimprovesthe productivityoftheop- erationcomparedtotheworst-casescenario.

At2.5h the conversionofcomponent Ais74%, andthe yield ofPis36%.Fig.12showstheestimateduncertainkineticparam- eters; basedon the update criterion (Eq. (19)andestimated val-

(10)

uesarewithintheworstcaseinterval)kineticparametersareﬁrst overwrittenat0.61h.Theaveragecomputationtimeis12.6sper iterations,hencethereal-timeoptimizationisfeasiblewiththisalgorithm.

Fig. 12 showshow the valuesof identified kinetic parameters gotothe realparameters asmoreinformationandmeasurement isavailable aboutthesystem. Lowrelative deviations(<1%)indi- catethattheidentifiedkineticparametersonlyslightlychange,so wecansaythattheidentifiedkineticparametersarenearthereal system,andwecanupdatetheuncertainparameters.

5. Conclusion

A framework to keep SBRs with exothermic reactions under controlinthewholeoperationusinganonlinearmodelpredictive control approach is proposed. The framework wastested on the semi-batch version of the Williams-Otto process including three reactions.Theproposed controlapproachcanalsohandletheun- certainkinetic parameters ofreactions.Inthiswork, theparame- tersofthefirstdominantreactionareconsideredthesourceofun- certaintyinthemodel.TheproposedframeworkconsistsofNMPC, EKF and an identification tool. The Modified Dynamic Condition was implemented into NMPC as an additional safety constraint, and the reactor temperature cannot exceed MAT. EKF is neces- sarytoestimatethestatevariablesofthereactorsystemandreac- tionrateconstants.Kineticparameterscanbeidentifiedwithleast squaresmethods basedon the estimatedreactionconstants after someformaltransformation.

We havecompared the Multi-Stage NMPC solved by the pro- gressivehedgingalgorithmandworst-casescenario.Eachresulted inaconservative solution,buttheworst-casescenarioNMPC has lowercomputationtime.InthecaseofMS-NMPC,thesize ofthe optimization problemincreases exponentially with the length of therobust horizonandwiththe uncertainparameters. Therefore, we have decided to extend the worst-case scenario NMPC with thestateestimationandidentiﬁcationalgorithms.Theresultsshow thattheproposedapproach canhandleuncertainkinetic parameters,andcanbe applied inreal reactorsystems inwhich reactor runawaycandeveloptoensuretheoptimalproduction.

Our future investigation will be about implementing the pro- posedcontrolapproachintoalab-scalereactorsystem.

DeclarationofCompetingInterest

Theauthorsdeclarethattheyhavenoknowncompetingﬁnan- cialinterestsorpersonalrelationshipsthatcouldhaveappearedto inﬂuencetheworkreportedinthispaper.

CRediTauthorshipcontributionstatement

Alex Kummer:Conceptualization, Methodology, Software,For- mal analysis, Investigation, Writing - original draft, Writing - review & editing, Visualization. Lajos Nagy: Conceptualization, Methodology, Writing - review & editing, Supervision. Tamás Varga:Conceptualization,Methodology,Writing-review&editing, Supervision.

Acknowledgement

We wouldliketoexpressouracknowledgementfortheﬁnan- cialsupportofSzéchenyi2020underGINOP-2.2.1-15-2017-00059. Tamás Varga’s contribution to this paper was supported by the UNKP-19-4 New National ExcellenceProgram of the Ministry for InnovationandTechnology.

References

Abel, O., Helbig, A., Marquardt, W., Zwick, H., Daszkowski, T., 20 0 0. Productivity optimization of an industrial semi-batch polymerization reactor under safety constraints. J. Process Control 10 (Aug. (4)), 351–362. doi: 10.1016/S0959-1524(99) 0 0 049-9 .

Albalawi, F., Durand, H., Christoﬁdes, P.D., 2018. Process operational safety via model predictive control: recent results and future research directions. Comput. Chem.

Eng. 114, 171–190. doi: 10.1016/j.compchemeng.2017.10.006 , Jun.

Arellano-Garcia, H., Barz, T., Dorneanu, B., Vassiliadis, V.S., 2020. Real-time feasi- bility of nonlinear model predictive control for semi-batch reactors subject to uncertainty and disturbances. Comput. Chem. Eng. 133, 106529. doi: 10.1016/j.

compchemeng.2019.106529 , Feb.

Ellis, M., Durand, H., Christoﬁdes, P.D., 2014. A tutorial review of economic model predictive control methods. J. Process Control 24 (Aug. (8)), 1156–1178. doi: 10.

1016/j.jprocont.2014.03.010 .

Fabiano, B., Vianello, C., Reverberi, A.P., Lunghi, E., Maschio, G., 2017. A perspec- tive on Seveso accident based on cause-consequences analysis by three different methods. J. Loss Prev. Process Ind. 49, 18–35. doi: 10.1016/j.jlp.2017.01.021 , Sep.

Findeisen, R. , Allgöwer, F. , Biegler, L.T.(Eds.) , 2007. Assessment and Future Directions of Nonlinear Model Predictive Control, 358. Springer Berlin Heidelberg, Berlin, Heidelberg .

Hall, R.J., “Explosion at T2 Laboratories,” 2010. https://www.aiche.org/academy/

webinars/explosion- t2- laboratories .

Holtorf, F., Mitsos, A., Biegler, L.T., 2019. Multistage NMPC with on-line generated scenario trees: application to a semi-batch polymerization process. J. Process Control 80, 167–179. doi: 10.1016/j.jprocont.2019.05.007 , Aug.

Kähm, W., Vassiliadis, V.S., 2018a. Thermal stability criterion integrated in model predictive control for batch reactors. Chem. Eng. Sci. 188, 192–207. doi: 10.1016/

j.ces.2018.05.032 , Oct.

Kähm, W., Vassiliadis, V.S., 2018b. Optimal Lyapunov exponent parameters for stability analysis of batch reactors with Model Predictive Control. Comput. Chem.

Eng. 119, 270–292. doi: 10.1016/j.compchemeng.2018.08.038 , Nov.

Kähm, W., Vassiliadis, V.S., 2019. Thermal stability criterion of complex reactions for batch processes. Chem. Eng. Res. Des. 150, 187–205. doi: 10.1016/j.cherd.2019.07.

028 , Oct.

Kühl, P. , Diehl, M. , Milewska, A. , Molga, E. , Bock, H.G. , 2007. Robust NMPC for a benchmark fed-batch reactor with runaway conditions. In: Findeisen, R., All- göwer, F., Biegler, L.T. (Eds.). In: Assessment and Future Directions of Nonlin- ear Model Predictive Control, 358. Springer Berlin Heidelberg, Berlin, Heidel- berg, pp. 455–464 .

Kanavalau, A., Masters, R., Kähm, W., Vassiliadis, V.S., 2019. Robust thermal stability for batch process intensiﬁcation with model predictive control. Comput. Chem.

Eng. 130, 106574. doi: 10.1016/j.compchemeng.2019.106574 , Nov.

Kao, C.-S., Hu, K.-H., 2002. Acrylic reactor runaway and explosion accident analysis. J. Loss Prev. Process Ind. 15 (May (3)), 213–222. doi: 10.1016/S0950-4230(01) 0 0 070-5 .

Khodadadi, H., Jazayeri-Rad, H., 2011. Applying a dual extended Kalman ﬁlter for the nonlinear state and parameter estimations of a continuous stirred tank reactor. Comput. Chem. Eng. 35 (Nov. (11)), 2426–2436. doi: 10.1016/j.compchemeng.

2010.12.010 .

Kummer, A., Varga, T., 2017. Feeding trajectory optimization in fed-batch reactor with highly exothermic reactions. Comput. Chem. Eng. 98, 1–11. doi: 10.1016/j.

compchemeng.2016.12.008 , Mar.

Kummer, A., Varga, T., 2019. Completion of thermal runaway criteria: two new criteria to deﬁne runaway limits. Chem. Eng. Sci. 196, 277–290. doi: 10.1016/j.ces.

2018.11.008 , Mar.

Kummer, A., Varga, T., Abonyi, J., 2019. Genetic programming-based development of thermal runaway criteria. Comput. Chem. Eng., 106582 doi: 10.1016/j.

compchemeng.2019.106582 , Sep.

Kummer, A., Varga, T., Nagy, L., 2020. Semi-batch reactor control with NMPC avoiding thermal runaway. Comput. Chem. Eng. 134, 106694. doi: 10.1016/j.

compchemeng.2019.106694 , Mar.

Lucia, S., Finkler, T., Engell, S., 2013. Multi-stage nonlinear model predictive control applied to a semi-batch polymerization reactor under uncertainty. J. Process Control 23 (Oct. (9)), 1306–1319. doi: 10.1016/j.jprocont.2013.08.008 .

Lucia, S., “Robust multi-stage nonlinear modelpredictive control,” 2014.

Mayne, D.Q., Seron, M.M., Rakovi ´c, S.V., 2005. Robust model predictive control of constrained linear systems with bounded disturbances. Automatica 41 (Feb. (2)), 219–224. doi: 10.1016/j.automatica.2004.08.019 .

Puschke, J., Mitsos, A., 2018. Robust feasible control based on multi-stage eN- MPC considering worst-case scenarios. J. Process Control 69, 8–15. doi: 10.1016/

j.jprocont.2018.07.004 , Sep.

Qin, S.J., Badgwell, T.A., 2003. A survey of industrial model predictive control tech- nology. Control Eng. Pract. 11 (Jul. (7)), 733–764. doi: 10.1016/S0967-0661(02) 00186-7 .

Qu, C.C., Hahn, J., 2009. Process monitoring and parameter estimation via unscented Kalman ﬁltering. J. Loss Prev. Process Ind. 22 (Nov. (6)), 703–709. doi: 10.1016/j.

jlp.2008.07.012 .

Rawlings, J.B. , Mayne, D.Q. , 2009. Model Predictive Control: Theory and Design. Wis:

Nob Hill Pub, Madison .

Rossi, F., Copelli, S., Colombo, A., Pirola, C., Manenti, F., 2015. Online model-based optimization and control for the combined optimal operation and runaway prediction and prevention in (fed-)batch systems. Chem. Eng. Sci. 138, 760–771.

doi: 10.1016/j.ces.2015.09.006 , Dec.

(11)

Rossi, F., Reklaitis, G., Manenti, F., Buzzi-Ferraris, G., 2016. Multi-scenario robust online optimization and control of fed-batch systems via dynamic model-based scenario selection. AIChE J. 62 (Sep. (9)), 3264–3284. doi: 10.1002/aic.15346 . Rossi, F., Manenti, F., Pirola, C., Mujtaba, I., 2016. A robust sustainable optimization &

control strategy (RSOCS) for (fed-)batch processes towards the low-cost reduction of utilities consumption. J. Clean. Prod. 111, 181–192. doi: 10.1016/j.jclepro.

2015.06.098 , Jan.

Rossi, F., Manenti, F., Buzzi-Ferraris, G., Reklaitis, G., 2017. Combined dynamic optimization, optimal control and online runaway detection prevention under uncertainty. Chem. Eng. Trans. 57, 973–978. doi: 10.3303/CET1757163 , 0.

Schweitzer, J.-M., López-García, C., Ferré, D., 2010. Thermal runaway analysis of a three-phase reactor for LCO hydrotreatment. Chem. Eng. Sci. 65 (Jan. (1)), 313–

321. doi: 10.1016/j.ces.2009.07.012 .

Scokaert, P.O.M., Mayne, D.Q., 1998. Min-max feedback model predictive control for constrained linear systems. IEEE Trans. Automat. Control 43 (Aug. (8)), 1136–

1142. doi: 10.1109/9.704989 .

Seki, H., Ogawa, M., Ooyama, S., Akamatsu, K., Ohshima, M., Yang, W., 2001. In- dustrial application of a nonlinear model predictive control to polymerization reactors. Control Eng. Pract. 9 (Aug. (8)), 819–828. doi: 10.1016/S0967-0661(01) 0 0 046-6 .

Sriram, M., Stevens, W.F., 1973. An example of the application of nonlinear programming to chemical-process optimization. Oper. Res. 21 (1), 296–304. doi: 10.1287/

opre.21.1.296 , Feb.

Thangavel, S., Lucia, S., Paulen, R., Engell, S., 2018. Dual robust nonlinear model predictive control: a multi-stage approach. J. Process Control 72, 39–51. doi: 10.

1016/j.jprocont.2018.10.003 , Dec.

van Woezik, B.A .A ., Westerterp, K.R., 20 0 0. The nitric acid oxidation of 2-octanol. A model reaction for multiple heterogeneous liquid–liquid reactions. Chem. Eng.

Process. 39 (Nov. (6)), 521–537. doi: 10.1016/S0255-2701(0 0)0 0 099-4 . Vanderbei, R.J., Shanno, D.F., 1999. An interior-point algorithm for nonconvex non-

linear programming. Comput. Optim. Appl. 13 (1/3), 231–252. doi: 10.1023/A:

1008677427361 .

Westerterp, K.R., Molga, E.J., 2004. No more runaways in ﬁne chemical reactors. Ind.

Eng. Chem. Res. 43 (Aug. (16)), 4585–4594. doi: 10.1021/ie030725m .

Yu, Z.J., Biegler, L.T., 2019. Advanced-step multistage nonlinear model predictive control: robustness and stability. J. Process Control 84, 192–206. doi: 10.1016/j.

jprocont.2019.10.006 , Dec.

Zhang, Z., Wu, Z., Durand, H., Albalawi, F., Christoﬁdes, P.D., 2018. On integration of feedback control and safety systems: analyzing two chemical process applications. Chem. Eng. Res. Des. 132, 616–626. doi: 10.1016/j.cherd.2018.02.009 , Apr.

Zhang, B., Hou, H., Hao, L., Zhu, J., Sun, Y., Wei, H., 2019. Identiﬁcation and optimization of thermally safe operating conditions for single kinetically controlled reactions with arbitrary orders in isoperibolic liquid-liquid semibatch reactors.

Chem. Eng. J. 375, 121975. doi: 10.1016/j.cej.2019.121975 , Nov.