Lyapunov exponents with Model Predictive Control for exothermic batch reactors

Proceedings, 10th IFAC International Symposium on Proceedings, 10th IFAC International Symposium on Advanced Control Chemical Processes Proceedings, 10th of IFAC International Symposium on at www.sciencedirect.com Available online Advanced Control Chemical Processes Proceedings, 10th of IFAC International Symposium on Shenyang,Control Liaoning, China, July 25-27, 2018 Advanced of Chemical Processes Shenyang, Liaoning, China, July 25-27, 2018 Advanced of China, Chemical Processes Shenyang,Control Liaoning, July 25-27, 2018 Shenyang, Liaoning, China, July 25-27, 2018

ScienceDirect

IFAC PapersOnLine 51-18 (2018) 423–428

Lyapunov exponents with Model Predictive Lyapunov Lyapunov exponents exponents with with Model Model Predictive Predictive Control for exothermic batch Lyapunov exponents with Modelreactors Predictive Control for exothermic batch reactors Control for exothermic batch reactors Control forahm, exothermic batch reactors Walter K¨ Dr Vassilios V. Vassiliadis *

Walter K¨ ahm, Dr Vassilios V. Vassiliadis * Walter K¨ ahm, Dr Vassilios V. Vassiliadis * Walter K¨ ahm, Dr Vassilios V. VassiliadisUniversity * Department of Chemical Engineering and Biotechnology, of Department of Chemical Engineering and Biotechnology, University of Department of Chemical Engineering and Biotechnology, University of Cambridge, Philippa Fawcett Drive, Cambridge, CB3 0AS, UK Cambridge, Philippa Fawcett Drive, Cambridge, CB3 0AS, UK Department of Chemical Engineering and Biotechnology, Cambridge, Philippa Fawcett Drive, Cambridge, CB3 University 0AS, UK of (e-mail: [email protected]). (e-mail: [email protected]). Cambridge, Philippa Fawcett Drive, Cambridge, CB3 0AS, UK (e-mail: [email protected]). * (e-mail: [email protected]). * (e-mail: [email protected]). (e-mail: [email protected]). * (e-mail: [email protected]). * (e-mail: [email protected]). Abstract: Abstract: Thermal Thermal runaways runaways cause cause significant significant safety safety issues issues and and financial financial loss loss for for industrial industrial Abstract: Thermal runaways cause significant safety issues andintensification financial lossofforprocesses industrial batch reactors due to the disruption of normal operation. The is batch reactors due to the disruption of normal operation. The intensification of processes is Abstract: Thermal cause significant safety issues andintensification financial lossoffor industrial batch reactors due torunaways the disruption ofcapable normal ofoperation. The processes is restricted, since control systems are not detecting stability boundaries of the system restricted, sincedue control systems are not capable ofoperation. detecting The stability boundariesofofprocesses the system batch reactors to the disruption of normal intensification is restricted, since control systems are For not this capable of detecting stability boundaries of the system and hence are overly conservative. purpose Lyapunov exponents are introduced as and hence since are overly conservative. For this purpose Lyapunov exponents are introduced as aa restricted, control systems are not capable of detecting stability boundaries of the system and hence are overly conservative. For this purpose Lyapunov exponents are introduced as of a stability criterion. It shown Lyapunov exponents can predict the stability criterion. It is is shown that that For Lyapunov exponents can correctly correctly predict the stability stability and hence aresystems. overly conservative. this purpose Lyapunov exponents are introduced as of a stability criterion. It is shown that Lyapunov exponents can correctly predict the stability of batch reactor This stability criterion is embedded in Model Predictive Control, which batch reactor systems. This stability criterion isexponents embeddedcan in Model Predictive Control, which stability criterion. It is shown that Lyapunov correctly predict the stability of batch reactor systems. This stability criterion is embedded in Model Predictive Control, which results in control scheme. This scheme allows controlled increase the results in a a novel novel control scheme. This scheme allows the the in controlled increase of of the reaction reaction batch systems. This stability criterion embedded Modeltime Predictive Control, which resultsreactor in a novel control scheme. This scheme allows the controlled increase of the reaction temperature to achieve a target conversion in aais reduced completion of reaction. temperature to achieve a target conversion in reduced completion time of reaction. results in a novel control scheme. This scheme allows the controlled increase of the reaction temperature to achieve a target conversion in a reduced completion time of reaction. © 2018, IFAC to (International Federation of Automatic Control) Hosting by Elsevier Ltd. All rights reserved. temperature achieve a target conversion in a reduced completion time of reaction. Keywords: Keywords: Process Process and and Control Control Monitoring, Monitoring, Batch Batch Process Process Modelling Modelling and and Control, Control, Process Process Keywords: Process and Control Monitoring, Batch Process Modelling and Control, Process Intensification Intensification Keywords: Process and Control Monitoring, Batch Process Modelling and Control, Process Intensification Intensification 1. INTRODUCTION Lyapunov exponents exponents (Strozzi (Strozzi and and Zaldvar, 1994) 1994) were later later 1. Lyapunov 1. INTRODUCTION INTRODUCTION Lyapunov exponents (Strozzi and Zaldvar, Zaldvar, 1994) were were later introduced to measure the stability of nonlinear systems introduced to the stability of nonlinear systems 1. INTRODUCTION Lyapunov exponents and Zaldvar, 1994) were later introduced to measure measure the theory. stability nonlinear systems using concepts from (Strozzi chaos Asofwill will be shown, shown, this using concepts from chaos theory. As be introduced to measure the stability of nonlinear systems using concepts from chaos theory. Asfor will be steady-state shown, this this criterion gives more reliable results non criterion gives more reliable results for non steady-state using concepts from chaos theory. Asfor will becost. shown, this criterion givescan more reliable results non steady-state systems, but impose high computational Batch processes processes are are very very important important for for the the chemical chemical indusindus- systems, but can impose high computational cost. criterion gives more reliable results for non steady-state Batch but can impose high computational cost. Batch processes are veryofimportant the chemical indus- systems, try, since since a vast vast amount amount specialityfor chemicals is produced produced In batch processes the of resystems, can impose highgeneration computational cost. try, a of speciality chemicals is In batch but processes the heat heat generation of exothermic exothermic reBatch processes are very important for the chemical industry, since a vast amount of speciality chemicals is produced in this reactor type. One advantage of using batch or In batch processes the heat generation of exothermic reactions decreases as the reagents are consumed. A timein this reactor type. One advantage of using batch or actions decreases as the reagents are consumed. A timetry,this sincereactor a vast amount of speciality chemicals isreactors, produced In batchdecreases processes heat generation of exothermic rein type.asOne advantage of using batch or semi-batch reactors, to is reagents are consumed. timedependent increaseasthe inthe reaction temperature in the theA stable stable semi-batch reactors, as opposed opposed to continuous continuous reactors, is actions increase in reaction temperature in in this reactor type. advantage of usingreactors, batch or actions decreases as the reagents are consumed. A timesemi-batch reactors, asOne opposed toand continuous is dependent the flexibility in process strategy control, due to the dependent increase in reaction temperature in the stable region leads leads to to an intensification intensification of of the process. process. The The ability ability the flexibility in process strategy control, due to semi-batch reactors, asan opposed toand continuous reactors, is region dependent increase inofreaction temperature inThe the imporstable the flexibility in process strategy and control, dueinvolves to the the possibility to control additional input. This region leadsinstability to an an intensification of the the isprocess. ability to detect the process of major possibility to control an additional input. This involves to detect instability of the process is of major importhe flexibility in process strategy and control, due to the region leads to an intensification of the process. The ability possibility to control an additional input. This involves processes operating operating in in non-steady non-steady state state operation, operation, which which to detect instability the systems. process is of major importance for the control of ofofsuch Furthermore, a reliprocesses for control systems. Furthermore, a possibility to control. controlin an additionalstate input. This involves to detect instability the process is of major imporprocesses operating non-steady operation, which tance are harder to tance for the the control ofofsuch such systems. Furthermore, a relireliable measure of stability can improve operator knowledge are harder to control. of can improve operator knowledge processes tancemeasure forand the therefore control of vastly such Furthermore, a reliare harderoperating to control.in non-steady state operation, which able able measure of stability stability can systems. improve operator knowledge greatly, reduce the risk of thermal Thermal runaways occur when when the the heat heat generated generated by by greatly, and therefore vastly reduce the risk of thermal are harderrunaways to control.occur able measure of stability can improve operator knowledge Thermal and therefore vastly reduce the risk of thermal runaways occurring in industry. This has an important Thermal runaways occur when the heatcapacity generated by greatly, runaways occurring in vastly industry. Thisthe hasrisk an of important exothermic reactions exceeds the cooling of the greatly, and therefore reduce thermal exothermic reactions exceeds the cooling capacity of the runaways occurring in industry. This has an important safety implication and decreases financial loss. Thermal runaways occur when the heat generated by exothermic reactionsresulting exceeds the cooling capacity of the safety reactor, ultimately ultimately in an an explosion. Stability implication and decreases financial loss. runaways occurring in industry. This has an important reactor, resulting in explosion. Stability safety implication and decreases financial loss. exothermic reactions exceeds cooling capacity ofisthe reactor, ultimately resulting in explosion. Stability criteria which reliably predict the howan close the system to Model Predictive Control Control (MPC)financial is an an advanced advanced control safety implication and decreases loss. criteria which reliably predict how close the system is Predictive (MPC) is control reactor, ultimately resulting in an explosion. Stability criteria which reliably predict how close the system is to to Model instability are therefore of major interest for the industry. Model Predictive Control (MPC) is an advanced control scheme which optimises the control variables of the system, instability are of interest for industry. which the control variables of criteria which reliably predict how close the system is to scheme Model considering Predictive Control (MPC) anInadvanced control instability are therefore therefore of major major interest for the the industry. scheme which optimises optimises theconstraints. control is variables of the the system, system, while system literature most Many criteria can been found in literature which describe while considering system constraints. In literature most instability are therefore of major interest for the industry. schemeconsidering which optimises theconstraints. control variables of the the system, Many criteria can been found in literature which describe while system In literature most MPC schemes implement a linearisation of system Many criteria can been found in literature which describe system stability. stability. The The first first stability stability criterion criterion was was developed MPC schemes implement aa linearisation of the system while considering constraints. In literature most system MPC schemes linearisation ofMPC the scheme system present, which implement cansystem be used used with a linear linear Many criteria canThe been found in literature which describe system stability. first stability criterionExplosion was developed developed by Semenov, who introduced the Thermal Thepresent, which can be with a MPC scheme MPC schemes implement a linearisation of the system by Semenov, who introduced the Thermal Explosion Thecan be2015). used with linear MPC scheme (Rawlingswhich and Mayne, Mayne, With asuch such formulation the system stability. The first stability criterion was by Semenov, who introduced the Thermal The- present, ory (Semenov, 1942). Introduction of many manyExplosion otherdeveloped stability (Rawlings and With aaa formulation the present, of which can be2015). used with asuch linear MPCtheoretischeme ory (Semenov, 1942). Introduction of other stability (Rawlings and Mayne, 2015). Withcan formulation the stability the closed-loop system be proven by Semenov, who introduced the Thermal Explosion Theory (Semenov, 1942). Introduction of many other stability criteria followed followed (Barkelew, (Barkelew, 1959) 1959) which, which, including including the the stability of the closed-loop system be theoreti(Rawlings Mayne, 2015). Withcan such formulation the criteria ofand the closed-loop system can bea proven proven theoretically by use of functions (DeHaan and Guay, ory (Semenov, 1942). Introduction ofwhich, manyofother stability criteria followed (Barkelew, 1959)bounds including the stability cally by the the use closed-loop of Lyapunov Lyapunovsystem functions (DeHaan and Guay, Thermal Explosion Theory, gave steady-state stability of the can be proven theoretiThermal Explosion Theory, gave bounds of steady-state cally by the use of Lyapunov functions (DeHaan and Guay, 2010). If no Lyapunov function can be found, end-point criteria followed (Barkelew, 1959) which, including the Thermal operating Explosion points. Theory, gave bounds of steady-state 2010). no Lyapunov function can be found, end-point cally byIf of Lyapunov functions (DeHaan and Guay, operating points. Ifthe nouse Lyapunov function can found,and end-point constraints are often employed. For be complex highly Thermal Explosion operating points. Theory, gave bounds of steady-state 2010). constraints are often employed. For complex and highly 2010). If no Lyapunov function can be found, end-point constraints are often employed. For complex and highly For dynamic systems one very common criterion is the nonlinear systems systems this this leads leads to to higher computational computational cost cost operating points. For dynamic one very criterion is constraints arehas often employed. For complex and frame. highly For dynamic systems systems one(Anagnost very common common is the the nonlinear nonlinear systems this leads to higher higher cost Routh-Hurwitz criterion and criterion Desoer, 1991). 1991). as the system to be simulated for aacomputational larger time as the system has to be simulated for larger time frame. Routh-Hurwitz criterion (Anagnost and Desoer, For dynamic one(Anagnost very common criterion is systhe as nonlinear this to higher computational cost Routh-Hurwitz criterion and Desoer, of 1991). the has to beleads simulated for acan larger time frame. This criterion systems only reliably describes the stability stability The usesystem ofsystems an online online stability criterion reduce the time This criterion only reliably describes the of sysThe use of an stability criterion can reduce the time Routh-Hurwitz criterion (Anagnost and Desoer, 1991). as the system has to be simulated for a larger time frame. This criterion only reliably describes the stability of sysof an criterion can system reduce the time tems at at steady-state. For For non steady-state steady-state systems systems this this The frameuse used byonline givingstability an indication of the stability tems used by giving an of stability This criterion onlyunreliable. reliably describes the stability of sysThe use of an stability criterion can system reduce the time tems at steady-state. steady-state. For non non steady-state systems this frame frame used byonline giving an indication indication of the the system stability criterion becomes at each point of the simulation. criterion becomes unreliable. each point the tems at steady-state. For non steady-state systems this at frame byof an indication of the system stability criterion becomes unreliable. at eachused point ofgiving the simulation. simulation. The divergence criterion quantifies how how much the the variables variables This work introduces stability criterion criterion based based on LyaLyacriterion becomes unreliable. at each point of the simulation. The divergence criterion quantifies work introduces aaa stability The divergence from criterion quantifies how much much the variables This This work introduces stability criterion basedit on on Lyaare ‘diverging’ each other based on the underlying punov exponents in a novel way, such that can be are ‘diverging’ from each other based on the underlying punov exponents in aaa novel way, such that can be The ‘diverging’ divergence criterion quantifies how much variables This work introduces stability criterion basedit are from each other based on thethe underlying exponents in Predictive novel way, such algorithms that it on canLyabea system equations. This definition comes from chaos theory punov integrated in Model Control in integrated in Model Predictive Control algorithms in a system equations. This definition comes from chaos theory are ‘diverging’ other based on thechaos underlying punov exponents in Predictive a studies novel way, such algorithms that caninbe system equations. definition comes from theory integrated in Model Control (Arnold, 1973), from butThis iseach not very reliable for non steady-state seamless manner. Case demonstrate the it efficacy ofa (Arnold, 1973), but is not very reliable for non steady-state seamless manner. Case studies demonstrate the efficacy of system equations. This definition theory seamless integratedmanner. in Model Control algorithms in ofa (Arnold, 1973), is not very forfrom non chaos steady-state CasePredictive studies demonstrate the efficacy systems and willbut therefore not reliable becomes considered further. systems will therefore not be (Arnold,and 1973), is not very for non further. steady-state seamless manner. Case studies demonstrate the efficacy of systems and willbut therefore not reliable be considered considered further. systems willIFAC therefore not beFederation considered further. Control) Hosting by Elsevier Ltd. All rights reserved. 2405-8963and © 2018, (International of Automatic

Copyright © 2018 IFAC 417 Copyright © under 2018 IFAC 417 Control. Peer review responsibility of International Federation of Automatic Copyright © 2018 IFAC 417 10.1016/j.ifacol.2018.09.337 Copyright © 2018 IFAC 417

2018 IFAC ADCHEM 424 Walter Kähm et al. / IFAC PapersOnLine 51-18 (2018) 423–428 Shenyang, Liaoning, China, July 25-27, 2018

the approach and the enhanced performance gained over more traditional PI control systems and MPC algorithms without such embedded stability criteria. 2. BATCH REACTOR SYSTEM In the model used for the subsequent simulations of batch processes, i.e. with constant volume, an irreversible, exothermic reaction is analysed which is given by: A+B → C (1) The model of the batch processes is based on differential equations for mass and energy balances. The reaction kinetics mainly consider component A, which are assumed to follow the Arrhenius equation (Davis and Davis, 2003). Examples of reactions with this kinetic scheme are polycondensation reactions, e.g. of dicarboxylic acid and diols, or the addition reaction for the synthesis of ethylene glycol from ethylene oxide and water. A diagram of the batch reactor system used is shown in Figure 1. Temperature transmitter

∆Hr is the enthalpy of the reaction, U is the heat transfer coefficient and A is the heat transfer area. A liquid, homogeneous reaction medium is modelled, which is why the physical properties and the volume of the reaction mixture are assumed to be constant. This leads to the following equation: V R ρ Cp

d (TR ) = r ([A] , TR ) (−∆Hr ) VR dt −U A (TR − TC )

(7)

The energy balance for the cooling jacket is given by: VC ρC CpC

d (TC ) = qC ρC CpC (TC,in − TC ) dt +U A (TR − TC )

(8)

where VC is the cooling jacket volume, ρC is the coolant density, CpC is the heat capacity of the coolant, TC,in and TC are the coolant inlet and cooling jacket temperature and qC is the volumetric coolant flow. The physical properties of the reaction mixture and of the reactor used in the simulations are shown in Table 1.

TT

Table 1. Process parameters for batch reactor simulations.

TIC

qC ; TC;in ; CpC ; ρC

qC ; TC

[A] ; [B] ; [C] V R ; T R ; Cp ; ρ

Fig. 1. Diagram of batch reactor with cooling jacket used for simulations. The mass balances for each reagent and product are: d [A] = −r ([A] , TR ) (2) dt d [B] = −r ([A] , TR ) (3) dt d [C] = +r ([A] , TR ) (4) dt   Ea n (5) r ([A] , TR ) = k0 [A] exp − R TR where [A], [B] and [C] are the concentrations of components A, B and C, t is time, r is the reaction rate, TR is the reactor temperature, k0 is the Arrhenius pre-exponential constant, n is the order of reaction with respect to component A, which is set to n = 2, Ea is the activation energy, and R is the universal molar gas constant. The energy balance for the reactor contents is given by: d (VR ρ Cp TR ) = r ([A] , TR ) (−∆Hr ) VR dt −U A (TR − TC ) (6) where VR is the reactor volume, ρ is the reactor content density, Cp is the heat capacity of the reactor contents, 418

Parameter

Value

VR

16 m3

Vc

1.2 m3

A

30.7 m2

TC, in

300 K

Ea /R

9525 K

k0

2.2 × 105 m3 kmol−1 s−1

qc, max

0.037 m3 s−1

ρ

1100 kg m−3

Cp

2330 J kg−1 K−1

ρC

1000 kg m−3

CpC

4180 J kg−1 K−1

U

600 W m−2 K−1

∆Hr

-75 ×106 J kmol−1

[A]0

13 kmol m−3

R

8.314 mol J−1 K−1

All simulations shown in this paper were carried out on a Dell XPS 13 with an Intel CoreTM i7-6560U processor, with operating system Windows 10 Home. The system dynamics were simulated using ode15s (Shampine et al., 1999) within MATLABTM . 3. LYAPUNOV EXPONENTS The Lyapunov exponents describe how state variables “drift off” after a large amount of time for an initial small perturbation, . The deviation of the state variables is assumed to follow an exponential profile, which enables to quantify if a stable system is present. A diagram showing the evolution of this deviation is shown in Figure 2.

2018 IFAC ADCHEM Walter Kähm et al. / IFAC PapersOnLine 51-18 (2018) 423–428 Shenyang, Liaoning, China, July 25-27, 2018

425

accuracy and numerical stability. If  was made smaller, the deviation inside the logarithm can get close to zero, therefore resulting in large negative numbers. The control variable is given by the governing equation (16) for PIcontrolled systems, and by 95% cooling for MPC controlled systems.

x(t) δx(t)

4. CONTROL SCHEMES

x(t0) + ǫ x(t0)

ǫ

4.1 PI Control structure

t0

t

The batch reactor temperature is controlled by varying the cooling water flow rate qC with a PI controller. The equation of the used PI controller is given by:

t=s

Fig. 2. Deviation of an initially perturbed state variable. In this case an unstable system is shown. The following expression quantifies the deviation of an initially perturbed state variable after time t:  exp (Λ (x0 ) t) =|x (t, x0 ) − x (t, x0 + ) | (9)   |x (t, x0 ) − x (t, x0 + ) | 1 (10) Λ (x0 ) = ln t  At the limit of a very small perturbation and infinite time:    δx (t, x0 )  1   (11) Λ (x0 ) = lim ln  t→∞ t δx0  This is known as the Lyapunov exponent (Strozzi and Zaldvar, 1994). Numerically, Lyapunov exponents can be evaluated by simulating several systems in parallel, for which each state variable is perturbed initially by an amount  = δx0 . Simulating the systems for an infinite amount of time is of course infeasible. Therefore a large time horizon is chosen instead, which is supposed to give a good approximation of the final value, known as the local Lyapunov exponent. This means that at each point in time, a long simulation has to be carried out in order to find the local Lyapunov exponent, given by:    δx (tf , x0 )  1  ln  Λl (x0 ) = (12)  tf δx0 The time tf in (12) is set in order to give a large time frame, which approximates the infinite horizon from the original definition. Other methods for evaluating the Lyapunov exponents are available (Melcher, 2003). Due to the heat generation and removal of the reaction the variables of interest are [A], TR and TC . Hence, the local Lyapunov exponents for these variables are evaluated by:    [A] (tf , [A]0 ) − [A] (tf , [A]0 + )  1  (13)  Λl,1 = ln   tf     TR (tf , TR,0 ) − TR (tf , TR,0 + )  1  (14)  Λl,2 = ln   tf     TC (tf , TC,0 ) − TC (tf , TC,0 + )  1  (15) Λl,3 = ln   tf 

The deviation of each variable is calculated at a final time of tf = 500 s, with an initial perturbation of  = 10−3 . Values for  in the range of 10−1 − 10−4 % of the variable analysed gave a good compromise between 419

u (t) = Kp (TR (t) − Tsp (t)) ˆt 1 (TR (t ) − Tsp (t )) dt + τI

(16)

t0

where Kp is the proportional parameter, τI is the integral parameter, Tsp (t) is the set point temperature at time t, u (t) is the control value and t is a dummy variable. The parameters of the PI controller are given in Table 2. The PI controller regulates the coolant flow, therefore controlling the reactor temperature. Table 2. Parameters for PI controller used in case studies. Parameter

Value

Proportional (P), Kp

10 m3 s−1 K−1

Integral (I), τI

1000 K s2 m−3

The current PI control is solely used to show how a stable system can become unstable due to a thermal runaway. The resulting system can be used to judge how well the Lyapunov exponents predict the system stability. No attempt was made to tune perfectly the PI controller since this is not the purpose of the following case study. 4.2 Model Predictive Control structure The analysis of stability for batch processes is incorporated into the classical MPC structure as shown in Figure 3. Predefined set-point

Stability analysis

Set- Optimal point control Controls algorithm

Process inputs

Process

Output

Fig. 3. Model Predictive Control scheme with integrated stability analysis. Model Predictive Control (MPC) is an advanced control scheme, in which an Optimal Control Problem is solved iteratively (Chuong La et al., 2017). The mathematical

2018 IFAC ADCHEM 426 Walter Kähm et al. / IFAC PapersOnLine 51-18 (2018) 423–428 Shenyang, Liaoning, China, July 25-27, 2018

formulation for MPC used in this work (Charitopoulos and Dua, 2016; Rawlings and Mayne, 2015) is given by: (17) min Φ (x (t) , u (t)) u(·)

subject to the system described in (2) − (8) and: ˆtf 2 Φ = (TR (t) − Tsp (t)) dt

uncontrolled rise in temperature, which accelerates the rate of reaction. Since an exothermic reaction is present, this leads to a thermal runaway. The temperature profile of this process is shown in Figure 5.

(18)

t0

Λl,i (tf ) ≤ 0 i = 1, 2, 3 (19) TR ≤ Tchem (20) (21) 0 ≤ qC ≤ qC,max t0 ≤ t ≤ tf (22) where t0 and tf are the initial time and final time of the simulation, and the chemical stability temperature is set to Tchem = 445 K. In the above problem formulation the coolant flow rate qC is equal to the control variable u (t), i.e. qC = u (t). The constraint in (19) ensures that the process does not enter an unstable region at the end of the horizon considered, i.e. at t = tf . This structure is different from that found in literature (Christofides et al., 2011), where Lyapunov functions are used. The problem given in (17) - (22) is solved using the SQP optimisation (Nocedal and Wright, 2006) algorithm within fmincon in MATLABTM . The algorithm proceeds with a “moving horizon”: At time t the optimal control action is evaluated for a control and prediction horizon of tcontrol and tprediction , respectively. This scheme is shown diagrammatically in Figure 4. Past control action Predicted control action Reference trajectory Measured process output Predicted model response

Fig. 5. Temperature profile of an exothermic batch reactor system. A step change in the set point temperature at t = 150 s leads to an unstable process. It can be seen clearly from Figure 5 that at time t ≥ 200 s the process turns unstable. A reliable stability criterion is required to identify this point of instability. For each 10 s interval the Lyapunov exponent values are evaluated. For a time frame of 500 s this gives 50 evaluations per simulation. The reliability to detect instability, as well as the computational time required to calculate each respective criterion are tested. The profiles for each Lyapunov exponent for the temperature profile given in Figure 5 are shown in Figure 6.

t=s tcontrol tprediction

Fig. 4. Diagram of Model Predictive Control with control and prediction horizons tcontrol and tprediction . The control action found by the optimisation algorithm is implemented only for the first step. After every iteration the algorithm is fed with new process data which, together with the included process model, lead to new predictions of the system behaviour. With this information the optimisation is carried out to find the optimal control values. 5. CASE STUDIES OF BATCH REACTIONS 5.1 PI Control An initially stable system is controlled using the PI control for the cooling jacket. At t = 150 s a step increase in the set point temperature of the reactor is implemented. This increase in the reactor temperature leads to an 420

Fig. 6. Lyapunov exponent profiles for a batch process going out of control. As can be seen in Figure 6 only the Lyapunov exponents predict the instability correctly. Before the new set point of 405 K is reached instability of the process is predicted. This instability is predicted at approximately 100 s. Therefore, the Lyapunov exponents are a reliable way of predicting instability of the batch process considered. One issue with Lyapunov exponents is the computational cost: per iteration the system model has to be simulated with initial perturbations in order to quantify stability. As the problem size increases, this can become an issue for online control applications. For this system, which evaluates the stability for three state variables, approximately 140 ms were required per iteration.

2018 IFAC ADCHEM Walter Kähm et al. / IFAC PapersOnLine 51-18 (2018) 423–428 Shenyang, Liaoning, China, July 25-27, 2018

5.2 Model Predictive Control Process intensification, which can be achieved by using a stability constraint with MPC, is demonstrated in this section. As a comparison a process controlled by MPC with Lyapunov exponents and MPC with a constant set point temperature are shown. The target conversion of the batch process is set to 85%. The maximum temperature allowed is set to Tchem = 445 K, which is the chemical stability of the process. The temperature profiles for each process are shown in Figure 7.

Fig. 7. Temperature profiles of batch processes controlled with MPC. For the process with constant set point of 405 K the temeprature increases beyond 550 K due to a thermal runaway. The conversion profiles for each process are shown in Figure 8.

427

To further illustrate the advantages of using Lyapunov exponents implemented with MPC, two different MPC strategies are presented below: • MPC with Lyapunov exponents, prediction horizon of tcontrol = 20 s, 1 step prediction • MPC without stability constraints, prediction horizon of tcontrol = 160 s, 8 step prediction with 20 s per step The performance for computational cost and stability are tested for the same batch reactor system, subject to an increase in set point temperature. The temperature profiles are shown in Figure 9.

Fig. 9. Temperature profiles of batch processes controlled by Model Predictive Control with and without Lyapunov exponent constraints. It can be seen from Figure 9 that the temperature profile for the MPC structure without stability constraints, but a larger control horizon, does not achieve the new set point temperature of 405 K in a stable manner. The resulting control profiles for the processes are shown in Figure 10.

Fig. 8. Conversion profiles of batch processes controlled with MPC. From Figure 8 it can be seen that the process including Lyapunov exponents as a stability constraint reaches the target conversion of 85% after 5,000 s in a controlled manner. The process with a constant temperature of 400 K needs 25,000 s. An increase of 5 K of the initial temperature already leads to an uncontrollable process. Therefore a significant decrease in reaction time (fivefold), whilst maintaining stability, was achieved with the inclusion of Lyapunov exponents as a measure of stability. This reduction of processing time results in earlier release of process units to be used for carrying out other tasks. On average, each iteration for the MPC scheme with Lyapunov exponents required 1.48 s. A control horizon of tcontrol = 20 s was used for both processes. The complete nonlinear process model was used for the control scheme, as the linearisation of this system can potentially lead to wrong predictions of the system behaviour. 421

Fig. 10. Temperature profiles of batch processes controlled with MPC. In Figure 10 it can be seen that the coolant flow for the process controlled by MPC with Lyapunov exponents increases rapidly just before t = 200 s as the boundary of instability has been reached. The Model Predictive Control structure including Lyapunov exponents as a constraint achieves a stable process, even though the control horizon is only 1/8th of the standard MPC formulation. The average computational cost for each MPC structure without the use of stability constraints and a control horizon of tcontrol = 160 s was 2.3 CPU seconds, whereas the MPC implementation using Lyapunov exponent constraints and a control horizon of tcontrol = 20 s required

2018 IFAC ADCHEM 428 Walter Kähm et al. / IFAC PapersOnLine 51-18 (2018) 423–428 Shenyang, Liaoning, China, July 25-27, 2018

an average of 2.0 CPU seconds. An MPC scheme with a control horizon of 500 s could be used, but with a larger horizon and more control intervals the optimisation problem increases, as well as the computational cost. Hence, the MPC structure using Lyapunov exponents as stability constraints not only improves the stability of the system, but also reduces the prediction horizon necessary and hence the computational cost of the optimisation algorithm. These are clear advantages to standard MPC strategies. This can lead to improved safety of operation for systems controlled by MPC in industry, which reduces financial loss due to interruptions of normal operation.

benefits over traditional control approaches, as well as the enhanced ability to intensify the underlying processes so as to achieve greater productivity. The implementation of such MPC schemes requires the direct use of highly nonlinear mechanistic process models for dynamic, i.e. non steady-state, systems because linearisation is not representative for fast dynamic regions. Hence many control steps are necessary for equivalent standard MPC schemes. Future work will continue by considering other, more complex reaction systems, as well as the impact of parametric uncertainty in the underlying process models used within the MPC scheme.

Problems still arise when using Lyapunov exponents as a stability criterion for nonlinear systems. The horizon over which the Lyapunov exponents are evaluated needs tuning, leading to varying results depending on the system. For exothermic batch reactions considered in this work Lyapunov exponents predicted the system stability reliably.

We thank the Engineering and Physical Sciences Research Council (EPSRC) and the Department of Chemical Engineering and Biotechnology, University of Cambridge, for funding the EPSRC PhD studentship for this project.

6. CONCLUSIONS AND FURTHER WORK

REFERENCES

A short review of common stability criteria for dynamic systems was presented. Key features of the Routh-Hurwitz criterion, Lyapunov exponents and the divergence criterion were given. The underlying theory of Lyapunov exponents was introduced and derived for the exothermic batch reactor system analysed. Additionally, the control systems used were explained and the implementation of Lyapunov exponents with Model Predictive Control (MPC), which is a novel control structure, was outlined. Advantages of using Lyapunov exponents as a measure of stability are the simple implementation and the reliability of the results obtained, once tuned correctly. Disadvantages include the need to tune the simulation horizon for Lyapunov exponents, as well as the computational cost as the problem size increases. From the MPC case studies it was shown that for unstable systems of small size, the MPC implementation using Lyapunov exponents resulted in a computationally cheaper and more stable process than the standard MPC implementation using a larger control horizon. More case studies of larger systems are required to prove this property. Furthermore, the use of Lyapunov exponents for a highly nonlinear system needs to be tested individually: there is no guarantee that Lyapunov exponents will give a correct prediction of system stability for every process. For batch reaction systems considered in this work Lyapunov exponents gave very reliable results. The use of Lyapunov exponents with MPC for an exothermic batch reaction lead to a profound decrease in reaction time. This is due to the capability of predicting the system stability along the process trajectory. Hence an intensification of the process was enabled while keeping the process under control at all times. This work has presented a totally new way of stabilising thermal runaway systems with an online MPC algorithm, while enhancing safety and performance of processes that can become unstable with detrimental effects leading to economic loss. The case studies presented demonstrate the 422

ACKNOWLEDGEMENTS

Anagnost, J. and Desoer, C. (1991). An elementary proof of the Routh-Hurwitz stability criterion. Circuits, Systems and Signal Processing, 10(1), 101–114. Arnold, V. (1973). Ordinary differential equations, chapter 3, 95–208. MA: MIT Press, Cambridge. Barkelew, C. (1959). Stability of chemical reactors. Chemical Engineering Progress Symposium Series, 25, 37–46. Charitopoulos, V.M. and Dua, V. (2016). Explicit model predictive control of hybrid systems and multiparametrix mixed integer polynomial programming. American Institute of Chemical Engineers Journal, 62, 3441–2460. Christofides, P.D., Liu, J., and Mu˜ noz de la Pe˜ na, D. (2011). Networked and Distributed Predictive Control, chapter 2, 13–45. Springer, London. Chuong La, H., Potschka, A., and Bock, H.G. (2017). Partial stability for nonlinear model predictive control. Automatica, 78, 14–19. Davis, M. and Davis, R. (2003). Fundamentals of Chemical Reaction Engineering, chapter 2, 53–56. McGraw-Hill. DeHaan, D. and Guay, M. (2010). Model Predictive Control, chapter 2, 26–58. Sciyo. Melcher, A. (2003). Numerische Berechnung der Lyapunov-Exponenten bei gew¨ ohnlichen Differentialgleichungen. Ph.D. thesis, Universit¨at Karlsruhe, Fakult¨ at f¨ ur Mathematik. Nocedal, J. and Wright, S. (2006). Numerical Optimization, chapter 18, 526–572. Springer, New York. Rawlings, J. and Mayne, D. (2015). Model Predictive Control: Theory and Design, chapter 1, 1–60. Nob Hill Publishing. Semenov, N. (1942). Thermal theory of combustion and explosion. Progress of Physical Science (U.S.S.R.), 23(3), 251–292. Shampine, L., Reichelt, M., and Kierzenka, J. (1999). Solving index-1 DAEs in MATLAB and Simulink. SIAM Review, 41, 538–552. Strozzi, F. and Zaldvar, J. (1994). A general method for assessing the thermal stability of batch chemical reactors by sensitivity claculation based on Lyapunov Exponents. Chemical Engineering Science, 49(16), 2681– 2688.