Parallelizable Real-Time Algorithm for Integrated Experiment Design MPC

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

ScienceDirect

IFAC PapersOnLine 51-18 (2018) 518–523

Parallelizable Parallelizable Real-Time Real-Time Algorithm Algorithm for for Parallelizable Real-Time Algorithm for Integrated Experiment MPC Parallelizable Real-Time Design Algorithm for Design MPC Integrated Experiment Integrated Experiment Design MPC Integrated Experiment Design MPC Xuhui Feng ∗∗∗ Yuning Jiang ∗∗∗ Mario Eduardo Villanueva ∗∗∗

Xuhui Feng ∗∗ Yuning Jiang ∗∗ Mario∗Eduardo Villanueva ∗∗ Xuhui Jiang Mario Boris Houska Xuhui Feng Feng Yuning Yuning Boris Jiang Houska Mario∗∗Eduardo Eduardo Villanueva Villanueva ∗ ∗ ∗Eduardo Villanueva ∗ Xuhui Feng ∗ Yuning Boris Jiang Mario Houska Boris Houska ∗ ∗ Houska ∗ of InformationBoris Science and Technology, ShanghaiTech ∗ School ∗ School of Information Science and Technology, ShanghaiTech ∗ School of Information Information Science and Technology, Technology, ShanghaiTech ShanghaiTech University, Shanghai, China ([email protected]) School of Science and ∗ University, Shanghai, China ([email protected]) School of Information and Technology, ShanghaiTech University, Shanghai, Science China ([email protected]) ([email protected]) University, Shanghai, China University, Shanghai, China ([email protected]) Abstract: Abstract: This This paper paper proposes proposes a a parallelizable parallelizable real-time real-time algorithm algorithm for for integrated integrated experimentexperimentAbstract: This paper proposes a parallelizable real-time algorithm for integrated design model predictive control (MPC). Integrated experiment design MPC is if Abstract: This paper proposes a parallelizable real-time algorithm for integrated experimentdesign model predictive control (MPC). Integrated experiment design MPC is needed neededexperimentif aa system system Abstract: This paper proposes a parallelizable real-time algorithm for integrated experimentdesign model predictive control (MPC). Integrated design MPC is needed system is not at and to on order able design model predictive control reference (MPC). Integrated experiment design MPC isin needed ifto system is not observable observable at a a tracking tracking reference and needs needs experiment to be be excited excited on purpose purpose in order if toaa be be able design modelthe predictive control (MPC). Integrated experiment design MPC is needed ifto a be system is not observable at a tracking reference and needs to be excited on purpose in order to be able to estimate system’s states and parameters. The contribution of this paper is a real-time is not observable at a tracking reference and needs to be excited on purpose in order able to estimate the system’s states and parameters. The contribution of this paper is a real-time is not observable at a tracking reference and needs to be excited on of purpose in order toreal-time be able to estimate the system’s states and parameters. The contribution of this paper is a real-time MPC algorithm using two processors. On the first processor an extended Kalman filter (EKF) to estimate the system’s states and parameters. The contribution this paper is a MPC algorithm using two processors. On the first processor an extended Kalman filter (EKF) to estimate the system’s states and parameters. The contribution of this paper is a real-time MPC algorithm using two processors. On the first processor an extended Kalman filter (EKF) as well as certainty-equivalent MPC which can MPC using two processors. On the firstcontroller processorare animplemented, extended Kalman (EKF) as wellalgorithm as a a parametric parametric certainty-equivalent MPC controller are implemented, whichfilter can provide provide MPC using two processors. On the firstcontroller processorare animplemented, extended Kalman (EKF) as as parametric certainty-equivalent which can immediate at sampling On second optimal design as well wellalgorithm as a a feedback parametric certainty-equivalent MPC are implemented, whichfilter can provide provide immediate feedback at high high sampling rates. rates. MPC On the thecontroller second processor, processor, optimal experiment experiment design as well as a feedback parametric certainty-equivalent MPC controller arecertainty-equivalent implemented, whichMPC can provide immediate feedback at high sampling rates. On the second processor, optimal experiment design (OED) problems are solved in parallel in order to perturb the control immediate at high sampling rates. On the second processor, optimal experiment design (OED) problems are solved in parallel in order to perturb the certainty-equivalent MPC control immediate feedback at high in sampling rates. On to theperturb second processor, optimal experiment design (OED) problems are solved in parallel in order to perturb the certainty-equivalent MPC control loop improving the accuracy of the state estimator at a lower sampling rate. We show that (OED) problems are solved parallel in order the certainty-equivalent MPC control loop improving the accuracy of the state estimator at a lower sampling rate. We show that this this (OED) problems are solved in parallel in order to perturb the certainty-equivalent MPC control loop improving the accuracy of the state estimator at a lower sampling rate. We show that framework can optimal objectives. approach is loop improving the accuracy of tradeoffs the state between estimatorOED at a and lowercontrol sampling rate. WeThe show that this this framework can achieve achieve optimal tradeoffs between OED and control objectives. The approach is loop improving the accuracy of tradeoffs the state between estimator at athat lower sampling rate. WeThe show that this framework achieve optimal OED and control objectives. approach is applied to biochemical process in to the proposed controller can achieve framework can achieve optimal OED and control objectives. The approach is applied to a acan biochemical processtradeoffs in order orderbetween to illustrate illustrate that the proposed controller can achieve framework can achieve optimal tradeoffs between OEDthat and the control objectives. The can approach is applied to a biochemical process in order to illustrate that the proposed controller can achieve superior control performance when compared to certainty-equivalent MPC. applied to a biochemical process in order to illustrate proposed controller achieve superior control performance when compared to certainty-equivalent MPC. applied a biochemical process in compared order to illustrate that the proposed controller can achieve superiorto control performance when compared to certainty-equivalent certainty-equivalent MPC. superior control performance when to MPC. © 2018, IFAC (International Federation of Automatic Hosting by Elsevier Ltd. All rights reserved. superior control performance when compared to Control) certainty-equivalent MPC. Keywords: Keywords: Model Model Predictive Predictive Control, Control, Parallelizable Parallelizable Algorithms, Algorithms, Optimal Optimal Experiment Experiment Design. Design. Keywords: Keywords: Model Model Predictive Predictive Control, Control, Parallelizable Parallelizable Algorithms, Algorithms, Optimal Optimal Experiment Experiment Design. Design. Keywords: Model Predictive Control, Parallelizable Algorithms, Optimal Experiment Design. 1. INTRODUCTION INTRODUCTION input in in order order to to increase increase the the accuracy accuracy of of future future state state and and 1. input 1. input to the accuracy of state and parameter estimates. Houska (2017) suggested a 1. INTRODUCTION INTRODUCTION input in in order order to increase increase the et.al. accuracy of future future state and parameter estimates. Houska et.al. (2017) suggested a selfselfINTRODUCTION input in order tocontroller, increase the accuracy of future state and Model predictive predictive 1. control (MPC) is is an an advanced advanced control control parameter estimates. Houska et.al. (2017) suggested aa selfreflective MPC which proceeds by minimizing minimizing parameter estimates. Houska et.al. (2017) suggested selfModel control (MPC) reflective MPC controller, which proceeds by Model (MPC) is control parameter estimates. Houska et.al. (2017) suggested a selftechnique, which control is capable capable of dealing dealing with inequality inequality MPC controller, which proceeds by minimizing the sum of a nominal control performance term as well as Model predictive predictive control (MPC) is an an advanced advanced control reflective reflective MPC controller, which proceeds by minimizing technique, which is of with the sum ofMPC a nominal control performance term as well as Model predictive (MPC) is anet.al. advanced control technique, which is of with inequality reflective controller, which byinherent minimizing stateand constraints (Mayne (2000); Rawlsum control performance term as as an additional term measuring theproceeds expected loss technique, which control is capable capable of dealing dealing with inequality the additional sum of of aa nominal nominal control performance term as well wellloss as stateand control control constraints (Mayne et.al. (2000); Rawl- the an term measuring the expected inherent technique, which(2009)). is capable of dealing with inequality stateand control constraints (Mayne et.al. (2000); Rawlthe sum of a nominal control performance term as wellloss as ings and Mayne Certainty-equivalent MPC con- an additional term measuring the inherent of of the controller (Stengel (1994)). state-and andMayne control constraints (Mayne et.al. (2000); Rawlan optimality additional term measuring the expected expected inherent loss of optimality of the controller (Stengel (1994)). ings (2009)). Certainty-equivalent MPC constateand control constraints (Mayne et.al. (2000); Rawlings and Mayne (2009)). Certainty-equivalent MPC conan optimality additional of term measuring the expected inherent loss trollers work well in practice if a reasonably accurate of the controller (Stengel (1994)). ings and Mayne (2009)). Certainty-equivalent MPC conof optimality of the controller (Stengel (1994)). trollers work well in practice if a accurate The above reviewed reviewed approaches have in in (1994)). common that that they they ings and Mayne (2009)). Certainty-equivalent MPC con- The trollers well in if aa reasonably reasonably accurate of optimality of the approaches controller (Stengel above have common model is work available (Zhu (2009)). Here, a basic basic requirement requirement trollers work well in practice practice if reasonably accurate model is available (Zhu (2009)). Here, a The above reviewed approaches have in common that they introduce additional hyperstates for predicting the accuThe above reviewed approaches have in common that they trollers work well in practice if a reasonably accurate model is available (Zhu (2009)). Here, a basic requirement introduce additional hyperstates for predicting the accuis that the controlled system is observable such that state model available (Zhu (2009)). Here, a basic requirement is that isthe controlled system is observable such that state introduce The above reviewed approaches common that they additional hyperstates for predicting accuracy of future future state estimates. estimates. Inhave particular, if the thethe variance introduce additional hyperstates for in predicting the accumodel isthe available (Zhu (2009)). Here, a basic requirement is that controlled system is observable such that state and parameter estimates can be determined online, e.g., by racy of state In particular, if variance is that the controlled system is observable such that state and parameter estimates can be determined online, e.g., by introduce additional hyperstates for predicting the accuracy of future state estimates. In particular, if the variance of future state estimates is penalized by augmenting the racy of future state estimates. In particular, if the variance is that the controlled system is observable such that state and parameter estimates can be determined online, e.g., by using an extended extended Kalman filter (EKF) or or moving moving horizon future state estimates is penalized by augmenting the and parameter estimates can be determined online, horizon e.g., by of using an Kalman filter (EKF) racy of future state estimates. Inmatrix particular, if the variance of future state estimates is penalized by augmenting the system dynamics with an EKF, EKF, valued hyperstates of future state estimates is penalized by augmenting the and parameter estimates can be determined online, e.g., by using an extended Kalman filter (EKF) or moving horizon system dynamics with an matrix valued hyperstates estimator (Diehl et.al. (2009); Rao et.al. (2001)). using an extended Kalman filter (EKF) moving horizon system estimator (Diehl et.al. (2009); Rao et.al.or(2001)). of future state estimates is penalized by augmenting the with EKF, matrix valued hyperstates have todynamics be included inan the MPC optimization problem, system dynamics with an EKF, matrix valued hyperstates using an extended Kalman filter (EKF) or moving horizon estimator (Diehl et.al. (2009); Rao et.al. (2001)). have to be included in the MPC optimization problem, estimator (Diehl et.al. (2009); Rao et.al. (2001)). fails to have system dynamics within anthe EKF, matrix valued hyperstates For some systems this observability requirement to be included MPC optimization problem, increasing the difficulty to solve it (Telen et.al. (2016)). have to be included in the MPC optimization problem, estimator (Diehl et.al. (2009); Rao et.al. (2001)). For some systems this observability requirement fails to increasing the difficulty to solve it (Telen et.al. (2016)). For this requirement fails have to bein included in the MPC(2016)) problem, hold if further is In the difficulty to solve it (Telen (2016)). Recently, (Feng and Houska tailored realFor some some systems this observability observability requirement fails to to increasing increasing the difficulty to solve itoptimization (Telenaa et.al. et.al. (2016)). hold if no no systems further precaution precaution is taken. taken. In such such scenarios, scenarios, in (Feng and Houska (2016)) tailored realForexample, some systems this observability requirement fails to Recently, hold if no further precaution is taken. In such scenarios, increasing the difficulty to solve it (Telen et.al. (2016)). for if the system is not observable at steady-state Recently, in (Feng and Houska (2016)) a tailored realtime optimization algorithm has been suggested, which hold if no further precaution is taken. In such scenarios, Recently, in (Feng algorithm and Houska (2016)) a tailored realtime optimization has been suggested, which for example, if the system is not observable at steady-state hold if no further precaution is observable taken. In such for if system is at Recently, in exploit (Feng algorithm and particular Houska a tailored realor ifexample, significant control excitations are needed needed toscenarios, improve time optimization has been attempts to the structure of MPC MPCwhich with for if example, if the the system is not not observable at steady-state steady-state time optimization algorithm has(2016)) been suggested, suggested, which or significant control excitations are to improve attempts to exploit the particular structure of with for example, ifofthe system isstate not observable at steady-state or if significant control excitations are needed to improve time optimization algorithm has been suggested, which the accuracy the online and parameter estimates, attempts to exploit the particular structure of MPC with integrated experiment design objectives. This approach or if significant control excitations are needed to improve attempts to exploit the particular structure of MPC with the of the online state andare parameter estimates, objectives. This approach or ifaccuracy significant control excitations needed to improve integrated the of the online state parameter attempts toexperiment exploit the design particular structure of MPC with it advisable the MPC objective function experiment design objectives. This approach improves the computation time of self-reflective MPC. theis accuracy of to theaugment online state and parameter estimates, integrated experiment design objectives. This approach it isaccuracy advisable to augment the and MPC objectiveestimates, function integrated improves the computation time of self-reflective MPC. the accuracy of to theaugment online state and parameter estimates, it is advisable the MPC objective function integrated experiment design objectives. This approach with an additional optimal experiment design (OED) obimproves the computation time of self-reflective MPC. it is advisable to augment the MPC objective function improves the time of self-reflective MPC. with additional optimal experiment design (OED) obTherefore, thecomputation current paper paper proposes an algorithm, algorithm, which it is an advisable to augment the MPC objective function with an additional optimal experiment design (OED) obimproves the computation time of self-reflective MPC. jective (Fisher (1935); Pukelsheim (1993)). During the last last the current proposes an which with an additional optimal experiment design (OED) ob- Therefore, jective (Fisher (1935); Pukelsheim (1993)). During the Therefore, the current paper proposes an algorithm, which splits the augmented MPC problem into a nominal and Therefore, the current paper proposes an algorithm, which with an additional optimal experiment design (OED) objective (Fisher (1935); Pukelsheim (1993)). During the last decade there have have beenPukelsheim numerous (1993)). suggestions on the howlast to splits the augmented MPC problem into aa nominal and jective (Fisher (1935); During decade there been numerous suggestions on how to Therefore, the current paperSection proposes an algorithm, which splits the augmented MPC problem into nominal and an experiment design part. 22 starts with a review splits the augmented MPC problem into a nominal and jective (Fisher (1935); Pukelsheim (1993)). During the last decade there have been numerous suggestions on how to an experiment design part. Section starts with a review integrate OED criteria in MPC (see Yan et.al. (2005)). For decade there beeninnumerous on howFor to an integrate OEDhave criteria MPC (see suggestions Yan et.al. (2005)). splits the augmented MPC problem22 starts into a with nominal and experiment design aa review of existing formulations of Section integrated experiment design an existing experiment design part. part. Section starts with review decade there have been numerous suggestions onaugment howFor to of integrate OED criteria in MPC (see Yan et.al. (2005)). example, Hovd and Bitmead (2005) suggest to formulations of integrated experiment design integrate OED criteria in MPC (see Yan et.al. (2005)). For example, Hovd and Bitmead (2005) suggest to augment an experiment design part. Section 2 starts with a review of existing formulations of integrated experiment design MPC. The main contribution of this paper is presented of existing formulations of integrated experiment design integrate OED criteria in MPC (see Yan et.al. (2005)). For example, and (2005) suggest the systemHovd dynamics by aa Riccati Riccati differential equation MPC. The main contribution of this paper is presented example, Hovd and Bitmead Bitmead (2005) differential suggest to to augment augment the system dynamics by equation of existing formulations integrated experiment design The main contribution of paper presented in Section 3, where we of propose to use use twois processors, MPC. The 3, main contribution of this this paper isprocessors, presented example, Bitmead (2005) differential suggeststate to augment the dynamics by aa Riccati implementing anand EKF. The matrix valued of this this MPC. in Section where we propose to two the system systemHovd dynamics byThe Riccati differential equation implementing an EKF. matrix valued stateequation of MPC. The main contribution of this paper is presented in Section 3, where we propose to use two processors, which communicate with each other and which provide in Section 3, where with we propose to use processors, the system dynamics byThe a Riccati differential implementing an EKF. matrix valued state this which communicate each other andtwo which provide Kalman filter used to variance of future implementing an be EKF. matrixthe valued stateequation of this which Kalman filter can can be usedThe to predict predict the variance ofof future in Section where sampling we propose toOn use communicate with each other and which provide feedback at 3, different rates. thetwo firstprocessors, processor which communicate with eachrates. other and which provide implementing anwhich, EKF. The matrix valued stateof of this Kalman filter can be used to predict the variance future state estimates, in turn, can be penalized in the feedback at different sampling On the first processor Kalman filter can be used to predict the variance of future state estimates, which, in turn, can be penalized in the which communicate with each other and which provide feedback at different sampling rates. On the first processor a parameterized standard MPC controller is running with feedback at different sampling rates. On the first processor Kalman filter can be used to predict the variance of future state which, in be the MPC objective. Similar Similar strategies to integrate integrate OEDin criteparameterized standard MPC controller is running with state estimates, estimates, which, strategies in turn, turn, can can be penalized penalized incritethe aafeedback MPC objective. to OED at different sampling rates. On the first processor parameterized standard MPC controller is running with high-sampling rate, while, on the other processor OED a parameterized standard MPC controller is running with state estimates, which, in turn, can be penalized in the MPC objective. Similar strategies to integrate OED critehigh-sampling rate, while, on the other processor OED ria in objective. MPC can can be be foundstrategies in (Heirung (Heirung et.al. (2012, (2012, 2015)). MPC Similar to integrate OED2015)). crite- high-sampling ria in MPC found in et.al. ahigh-sampling parameterized standard MPC controller is running with rate, while, on the other processor OED problems are solved. An associated communication scheme rate, while, on the other processor OED MPC objective. Similar strategies to integrate OED criteria in MPC can be found in (Heirung et.al. (2012, 2015)). problems are solved. An associated communication scheme Notice that there is also a number of articles on perria in MPC be found (2012, on 2015)). Notice that can there is alsoina(Heirung number et.al. of articles per- problems high-sampling rate, while, on the other processor OED are solved. An associated communication scheme uses ideas from the field of augmented Lagrangian based problems are solved. An associated communication scheme ria in MPC can beMPC found inaa(Heirung et.al. (2012, on 2015)). Notice that there is number of articles persistently exciting (Hernandez and Trodden (2015); ideasare from the field of augmented Lagrangian based Notice that there is also also numberand of Trodden articles on per- uses sistently exciting MPC (Hernandez (2015); problems solved. An associated communication the of Lagrangian based alternating direction inexact Newton (ALADIN) methuses ideas ideas from from the field field of augmented augmented Lagrangianscheme based Notice that there is also a number of(2015); articlesZacekova on per- uses sistently exciting MPC (Hernandez and Trodden (2015); Marafioti et.al. (2014); Mesbah et.al. alternating direction inexact Newton (ALADIN) methsistently exciting MPC (Hernandez and Trodden (2015); Marafioti et.al. (2014); Mesbah et.al. (2015); Zacekova uses ideas from the(2016)), field of which augmented Lagrangian based alternating direction inexact Newton (ALADIN) methods (Houska et.al. ensures that both proalternating direction inexact Newton (ALADIN) methsistently exciting MPC (Hernandez and Trodden (2015); Marafioti et.al. (2014); Mesbah et.al. (2015); Zacekova ods (Houska et.al. (2016)), which ensures that both proet.al. (2013)), which haveMesbah appearedet.al. recently and Zacekova which all all alternating direction inexact Newton (ALADIN) methMarafioti et.al. (2014); (2015); et.al. (2013)), which have appeared recently and which ods (Houska et.al. (2016)), which ensures that both processors can find a compromise between control excitation ods (Houska et.al. (2016)), which ensures that both proMarafioti et.al. (2014); Mesbah et.al. (2015); Zacekova et.al. (2013)), which have appeared recently and which all cessors can find a compromise between control excitation propose variants to perturb the nominally optimal control et.al. (2013)), which have appeared recently and which all cessors propose variants to perturb the nominally optimal control ods (Houska et.al. (2016)), which ensures that excitation both procan aa compromise between control can find find compromise between control excitation et.al. (2013)), which have appeared recently and which all cessors propose variants to the optimal control propose variants to perturb perturb the nominally nominally optimal control cessors can find a compromise between control excitation propose to perturb the nominally optimal control 2405-8963variants © 2018, IFAC (International Federation of Automatic Control) Hosting by Elsevier Ltd. All rights reserved.

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

2018 IFAC ADCHEM Shenyang, Liaoning, China, July 25-27, 2018Xuhui Feng et al. / IFAC PapersOnLine 51-18 (2018) 518–523

519

for the purpose of improving future state estimates and nominal control performance. Section 4 compares the runtime and control performance of the algorithm for nominal and integrated-experiment design MPC by studying a challenging biochemical process that is not observable at its steady-state. Section 5 concludes the paper.

Remark 2. Notice that there exist many articles on how to formulate MPC problems of the above form such that suitable closed-loop stability or other closed-loop performance criteria are met (Rawlings and Mayne (2009)).

2. INTEGRATED OED BASED MPC

For many processes the above certainty equivalent MPC in combination with an EKF works reasonably well in practice and, of course, in this case, no further modifications are needed. However, in some situations, e.g., if the MPC controller attempts to track a steady state at which the system states are not observable, EKF-MPC cascades may fail to perform well (Houska et.al. (2017); Telen et.al. (2016)). In this case, the predicted variance of future state estimates can be penalized in the objective of the MPC controller. This leads to an augmented optimization problem of the form

This section reviews methods for including OED criteria in the objective of an MPC controller. The focus is on discrete-time control systems of the form xk+1 = f (xk , uk ) + wk ηk = h(xk ) + vk ,

(1)

where xk ∈ Rnx denotes the state, uk ∈ Rnu the control input, and ηk ∈ Rnh the measurement at time k. The functions f : Rnx × Rnu → Rnx and h : Rnx → Rnh are assumed to be nonlinear and at least three times continuously differentiable. We are interested in analyzing OED criteria since the process noise wk and the measurement error vk are random variables. The system states have to be estimated frequently from the most recent measurements. As the control uk enters f nonlinearly, the control input might influence the observability properties of the system. In the following, it is assumed that the vk s and wk s are uncorrelated in time with mean 0 and given variances h V ∈ Sn++ and W ∈ Sn+x , respectively 1 . Remark 1. Notice that in the above framework uncertain parameters can be included by introducing auxiliary state variables that are invariant with respect to k, i.e., by stacking equations of the form pk+1 = pk to the discrete time recursion and regarding these parameters as states. 2.1 Model Predictive Control

l(xk , uk ) + m(xN ) k=0    ∀k ∈ {0, . . . , N − 1}, xk+1 = f (xk , uk ), x0 = x ˆ, s.t.   c(xk , uk ) ≤ 0 .

L(u, x ˆ) = min x

min

x,u,Σ

s.t.

N −1 

l(xk , uk ) + m(xN ) +

N 

Φ(xk , uk , Σk )

k=0 k=0  ∀k ∈ {0, . . . , N − 1},      xk+1 = f (xk , uk ), x0 = x ˆ

(3)

ˆ  Σk+1 = F (xk , uk , Σk ), Σ0 = Σ     c(xk , uk ) ≤ 0 ,

where Σ denotes the matrix-valued state of the EKF with F (x, u, Σ) := A(x, u)G(x, u, Σ)A(x, u)T + W

Here,

G(x, u, Σ) :=  −1 Σ − ΣC(x)T C(x)ΣC(x)T + V C(x)Σ .

(4)

∂f ∂h (x, u) and C(x) = (x) ∂x ∂x denote the first order derivatives of f and h with respect to the states. A(x, u) =

Certainty-equivalent MPC proceeds by solving online optimization problems of the form N −1 

2.2 Integrated Optimal Experiment Design

(2)

Here, l : Rnx × Rnu → R denotes the stage cost and c : Rnx × Rnu → R the mixed control- and state constraints. The prediction is computed under a certaintyequivalence paradigm, i.e., as if the current state estimate x ˆ was known exactly and as if there was no process noise, wk = 0. Of course, in practice, the certainty-equivalent controller only achieves reasonable control performance if the state estimates, e.g., found by an EKF, are reasonably accurate and if the controller is run in a receding horizon way; that is, the above problem is re-solved whenever a new state measurement x ˆ is available and the corresponding first element of the control sequence, denoted by u∗0 , is sent to the real process. 1

The set of real symmetric and positive semidefinite n×n matrices is n denoted by Sn + . Similarly, S++ denotes the set of symmetric positive definite matrices.

513

Notice that there are many ways to choose the functions Φ : Rnx × Rnu × Sn+x → R. For example, in Hovd and Bitmead (2005) it is suggested to use a weighted trace Φ(xk , uk , Σk ) := Tr (αk Σk ) , where αk ∈ Sn+x can be any scaling matrices that the user can choose. In principle, one could also use other standard OED criteria replacing the trace with the maximum eigenvalue (E-criterion) or a determinant (D-criterion) (Ljung (1999)). In general, the function Φ could even depend on xk and uk . This is for example, the case if a selfreflective MPC controller is used (Houska et.al. (2017)), which chooses Φ in such a way that the additional term in the objective of (3) can be interpreted as the controller’s own expected loss of optimality in the presence of the random process noise wk and the random measurement errors vk . For an overview of different choices of Φ in the context of MPC we refer to Telen et.al. (2016). Notice that all these problem formulations can be written in the form (3) by defining the function Φ appropriately. In order to avoid misunderstanding, given the large amount of articles about integrated experiment design based MPC, problem formulation (3) is not an original contribution of this paper. Rather, the contribution of this

2018 IFAC ADCHEM 520 Shenyang, Liaoning, China, July 25-27, 2018Xuhui Feng et al. / IFAC PapersOnLine 51-18 (2018) 518–523

paper is the development of an algorithm that can solve problems of the form (3) efficiently and in real-time. Notice that such algorithms are of practical relevance, since the additional matrix-valued state Σ makes (3) much more difficult to solve than the original MPC problem (2). Thus, in order to make MPC with additional OED objectives useful for practical applications, tailored algorithms and software are needed. The goal of the following section is to propose such an algorithm that exploits the structure by alternating between MPC and OED objectives. 3. ALTERNATING EXPERIMENT DESIGN AND CONTROL OBJECTIVES

x

s.t.

N 

Φ(xk , uk , Σk )

k=0 

∀k ∈ {0, . . . , N − 1},      xk+1 = f (xk , uk ), x0 = x ˆ

Initialization: Initial guesses for z and λ, positive definite scaling matrices HL and HE , and a tuning parameter ρ > 0. Schedule of Processor 1: Step 1a (W): ˆ as sent by Processor 2. Wait for x ˆ and Σ

This section proposes a real-time algorithm for solving (3). Here, the main idea is to split the overall problem into two interconnected optimization problem that can be solved in parallel and on two processors. The first processor solves certainty-equivalent MPC problems at a high sampling rate while the second processor solves OED problems at a lower sampling rate, in order to support the decisions of the other processor. In order to explain how this works, we introduce the shorthand notation ˆ := min E(u, x ˆ, Σ)

Algorithm 1: Alternating OED-MPC

Step 1b (OED): Solve the augmented OED problem ˆ + λT uE + ρ uE − z2 . min E(uE , x ˆ, Σ) 2 uE 2 and set gE = ρ(z − uE ) − λ.

(7)

Step 1c (W): Wait for the next update of uL as sent by Processor 2 and set gL = ρ(z − uL ) + λ. Step 1d (c-QP):

(5)

ˆ  Σk+1 = F (xk , uk , Σk ), Σ0 = Σ     c(xk , uk ) ≤ 0 . Notice that the function E is more expensive to evaluate than L, as the additional matrix-valued state Σ needs to be propagated forward in time. Next, the original OED based MPC problem (3) can be written in the equivalent consensus form   ˆ ˆ) + E(uE , x ˆ, Σ) s.t. uL = uE . (6) min L(uL , x uL ,uE

The above formulation reveals a parallelizable structure, which suggests to use parallelizable control algorithms. Notice that there exists mature MPC software (see Houska et.al. (2011)), as well as OED algorithms and softwares, see, e.g., Bauer et.al. (2000). Thus, it is desirable to make use of these existing software packages by developing a splitting scheme that exploits the structure of (6). However, OED problems are (at least without further reformulation) typically non-convex (Pukelsheim (1993)). Consequently, convex splitting methods such as the alternating direction method of multipliers (ADMM) cannot be applied. Therefore, the following algorithmic developments are based on the ALADIN method (Houska et.al. (2016)), as outlined next. 3.1 ALADIN as an OED-MPC splitting algorithm By applying ALADIN to (6), a distributed algorithm is obtained that solves decoupled OED and MPC problems in parallel. Algorithm 1 presents the main algorithmic steps for a real-time variant of ALADIN, which, in contrast to the standard ALADIN algorithm, solves multiple instances of the MPC problem at a higher frequency 514

Solve the consensus QP 1 1   min ∆u L HL ∆uL + ∆uE HE ∆uE + gL ∆uL ∆uE ,∆uL 2 2  (8) +gE ∆uE + λ ∆uL s.t.

uL + ∆uL = uE − ∆uE

| λQP ,

set z = uL + ∆uL and λ = λQP , send (z, λ) to Processor 2, and go back to Step 1a. Schedule of Processor 2: Step 2a (EKF): Wait for the measurement η and update the EKF in order to compute a new state estimate x ˆ and an ˆ associated variance Σ. Step 2b (MPC): If an update for z and λ (from Processor 1) is available, update these variable. Next, solve ρ 2 min L(uL , x ˆ) − λT uL + uL − z2 , (9) uL 2 send the first element of the optimal control u0 = (uL )0 to the real process, send uL to Processor 1, and continue with Step 2a. in order to give fast feedback. The schedule of the two processors from Algorithm 1 is additionally visualized in Figure 1 highlighting that Processor 1 and Processor 2 are running their main steps in parallel. Problem (9) remains an MPC optimization problem in standard form, but augmented albeit with the state and terminal cost. Thus, it can be solved with standard MPC solvers, where the additional linear terms as well as the additional least squares term hardly introduce any additional numerical

2018 IFAC ADCHEM Shenyang, Liaoning, China, July 25-27, 2018Xuhui Feng et al. / IFAC PapersOnLine 51-18 (2018) 518–523

521

Fig. 1. Schedule for the main steps of Algorithm 1. Processor 1 and Processor 2 run in parallel. difficulties. Consequently, Processor 2 can run with the same sampling time as a standard MPC controller for the same system, which should be considered as one of the main advantages of Algorithm 1 compared to other methods, which solve problem (3) directly.

proposed parallelizable algorithm can also be in combination with the sigma-point approach (see Kawohl (2007)), which approximates Σ more accurately.

Algorithm 1 can be interpreted as a tailored variant of the generic ALADIN algorithm that has been proposed in (Houska et.al. (2016)). Consequently, the local contraction proof of the ALADIN iterates from Houska et.al. (2016) can be applied one-to-one, thus the iterates of Algorithm 1 contract in a local neighborhood of a primal dual solution (u∗ , λ∗ ) solution of (6), i.e., as long as (z, λ) ≈ (u∗ , λ∗ ), and under the assumption that there is no noise present (see also Feng and Houska (2016)). In this paper, the matrices HL and HE are considered as scaling matrices, which can in practice be pre-computed by using a pre-conditioner, and ρ is kept constant. More details regarding how to choose the scaling and tuning factors of ALADIN can be found in Houska et.al. (2016). Of course, in practice, Algorithm 1 has to be run in the presence of process noise, but an empirical observation is that Algorithm 1 still controls the system reasonably well in this case, although a more detailed stability and robustness analysis of this empirical observation is beyond the scope of this paper. Remark 3. As mentioned above, the run-time complexity of most solvers that can deal with problem (9) is of order O(N n3x ), as this problem has the same size as standard MPC problems. On the other hand, solving problem (7) with a generic solver usually leads to an algorithm with run-time complexity O(N n6x ) as the Kalman filter states are matrix-valued. However, it turns out that one can develop tailored algorithm that can solve this problem with run-time complexity O(N n3x ), too, by exploiting the particular structure of the algebraic Riccati recursion (Bittanti et.al. (1991)), as for example discussed by Telen et.al. (2013) or also, in another variant, in Feng and Houska (2016). If such advanced OED solvers are used, the OED solver takes in practice still longer than the MPC solver, but the run-time ratio between the two solvers does not depend on the number of states. Remark 4. One could also imagine real-time variants of Algorithm 1. E.g., instead of solving the decoupled augmented OED and MPC problems to optimality one could apply one SQP step per iteration following the classical real-time scheme as pioneered by Diehl et.al. (2002) in order to solve the decoupled NLPs only approximately. Remark 5. Notice that the EKF may be inaccurate for nonlinear system in the presence of large noise, but the

4.1 Dynamic system model

515

4. NUMERICAL RESULTS

We consider a controlled chemical reactor with given dynamics ∀t ∈ [0, T ], χ(t) ˙ = g(χ(t), ν(t)), χ(0) = x , (10) where its right-hand expression is given by   −(D + k1 )χ1 + k2 χ2 χ3 + ν1   g(χ, ν) =  −Dχ2 − k2 χ2 χ3 + k1 χ1 + ν2  .

−Dχ3 − k2 χ2 χ3 + ν3 Here, χ [ Lg ] denotes concentrations of three substances, ν g ] the feeding rates, k1 and k2 [ 1s ] the reaction rates, and [ L·s D [ 1s ] the dilution rate. In the following, f (x, u) = χ(τ ) denotes the solution of the associated ODE for a small step-size τ > 0 and a piecewise constant control input. In this case study, a Runge-Kutta integrator of order 4 is used to evaluate the function f with high numerical accuracy. Moreover, the stage cost l, the Mayer term m and the constraint c are given by 1 1 l(xk , uk ) = xk − xref 2Q + uk − uref 2R , 2 2 1 m(xN ) = xN − xref 2QN , 2 as well as c(xk , uk ) = u − uk . Here, we assume that only the first state can be measured, h(x) = x1 . The numerical values for the MPC horizon length, references, and other parameters can be found in Table 1. 4.2 Implementation details Algorithm 1 for this particular case study has been implemented by using the algorithmic differentiation software CasADi (Andersson et.al. (2012)) in combination with automatic C-code generation. Moreover, we use the software qpOASES as a QP solver (Ferreau et.al. (2008)). All the results below are obtained on a macOS Sierra operating system with two 3.3 GHz Intel Core i7 processors and 16 GB, 2133 MHz LPDDR3. 4.3 Control performance Figure 2 shows a comparison of the closed-loop trajectories that are obtained by running Algorithm 1 and certainty-

2018 IFAC ADCHEM 522 Shenyang, Liaoning, China, July 25-27, 2018Xuhui Feng et al. / IFAC PapersOnLine 51-18 (2018) 518–523

Fig. 2. [Discrete-time system] Closed-loop state trajectories x with random noise as well as the associated closedloop input profile u obtained by running the alternating OED-MPC Algorithm 1 (red, dotted), by running the certainty-equivalent MPC (blue, dotted), and the reference states xref and controls uref (black, solid). Name

 2 M 1  OED-MPC OED-MPC g l(xi , ui ) ≈ 0.53 M L2 k=1

Symbol

Value

dilution rate

D

0.1 [ 1s ]

reaction rates

k1 , k 2

0.1, 0.5 [ 1s ]

discrete-time step-size

τ

0.5 [s]

MPC horizon length

N

10

initial state estimate

yˆ

g [1.0, 5.0, 0.0]T [ L ]

initial state variance

 Σ

g diag (0, 0, 0) [ L 2]

xref

g [1.0, 5.0, 0.0]T [ L ]

control reference

uref

g [0.6, 0.0, 0.0]T [ L·s ]

lower control bound

u

g [0, −1, −1]T [ L·s ]

measurement error variance

V

g 2.5 · 10−4 [ L 2]

process noise variance

W

g diag (0, 2.56, 0) [ L 2]

state weighting matrix

Q, QN

diag (1, 1, 1)

penalty term

R

4.4 Run-time performance

2

state reference

control weighting matrix

for randomly generated uniformly distributed uncertainty. Consequently, the proposed integrated OED-MPC controller achieves over three times better performance than certainty-equivalent MPC.

Although the weighted A-criterion comprises 18 states and 3 controls in total for this particular case, the runtime is in the microsecond range by combining the proposed parallelizable algorithm with the acceleration scheme from Feng and Houska (2016), which exploits the particular structure of integrated experiment design problem.

2

diag (1, 1, 5)

ρ

2

Processor 1 (OED processor)

[s2 ]

103

Table 1. Parameter values. equivalent MPC (solely running the control solver without any assistance from the OED solver) on the above case study, respectively. Here, the function Φ is a weighted trace term using the self-reflective weighted A-criterion as proposed in Houska et.al. (2017). Notice that for this particular case study the system is not observable at its steady states xref , since only the first state component of the system can be measured. Consequently, the controller finds an optimal tradeoff between estimation accuracy and tracking performance by steering the system to a steady-state that is observable, but not exactly equal to xref , as one would expect from an integrated OED based MPC controller. Their average control performance illustrates the difference between two controllers. By performing the closed-loop simulation for a sufficiently long time M , the average performance of the certainty-equivalent MPC controller is  2 M g 1  MPC MPC l(xi , ui ) ≈ 1.62 M L2 k=1

and the corresponding value of performing Algorithm 1 is 516

CPU time [µs]

%

Step 1a (W)

−

−

Step 1b (OED)

40

42

Step 1c (W)

−

−

Step 1d (c-QP)

55

58

Total time

95

100

Processor 2 (control processor)

CPU time [µs]

%

Step 2a (EKF)

≤1

≈4

Step 2b (MPC)

24

96

Total time

25

100

Table 2. Run-time associated with the different steps of Algorithm 1 for both processors. Table 2 summarizes the run-time of different steps of Algorithm 1. The decoupled problems (7) and (9) are solved by using real-time iterations as discussed in Remark 4. The time for solving both the OED and the coupled QP is about 95µs, which corresponds to around four times of the sampling time of the certainty-equivalent MPC loop. The waiting times in Step 1a and 1c were set to 0. 5. CONCLUSION This paper has presented the alternating MPC-OED Algorithm 1, which finds optimal tradeoffs between OED

2018 IFAC ADCHEM Shenyang, Liaoning, China, July 25-27, 2018Xuhui Feng et al. / IFAC PapersOnLine 51-18 (2018) 518–523

criteria and MPC control performance objectives. The main contribution of this tailored algorithm is that it uses two processors: Processor 2 runs at the same sampling rate as standard MPC without additional OED objectives would do. Processor 1 is used to support the decisions of Processor 2 at a lower sampling rate by perturbing the gradient of the MPC controller in order to improve the accuracy of the state estimates in the presence of random process and measurement noise. A case study has illustrated the practical advantages of this algorithm in terms of both run-time as well as control performance. A more detailed robustness and stability analysis supporting the algorithmic developments will be part of future work. ACKNOWLEDGEMENTS The authors gratefully acknowledge the financial support by the National Science Foundation China (NSFC), Nr. 61473185, as well as ShanghaiTech University, Grant-Nr. F-0203-14-012. REFERENCES Andersson, J., Akesson, J., and Diehl, M. (2012). CasADi: A Symbolic Package for Automatic Differentiation and Optimal Control. Recent Advances in Algorithmic Differentiation, 297–307, Springer, 2012. Bauer, I., Bock, H.G., K¨ orkel, S., and Schl¨ oder, J.P. (2000). Numerical methods for optimum experimental design in DAE systems. Journal of Computational and Applied Mathematics, 120: 1–25, 2000. Bittanti, S., Colaneri, P., and De Nicolao, G. (1991). The periodic Riccati equation. In Willems Bittanti, Laub, editor, The Riccati Equation. Springer, 1991. Diehl, M., Bock, H.G., Schl¨ oder, J.P., Findeisen, R., Nagy, Z., and Allg¨ ower, F. (2002). Real-time optimization and Nonlinear Model Predictive Control of Processes governed by differential-algebraic equations. Journal of Process Control, 12(4): 577–585, 2002. Diehl, M., Ferreau, H.J., and Haverbeke, N. (2009). Efficient Numerical Methods for Nonlinear MPC and Moving Horizon Estimation, Nonlinear Model Predictive Control, 384: 391–417, 2009. Feng, X. and Houska, B. (2016). Real-time algorithm for self-reflective model predictive control. Journal of Process Control, 2017, https://doi.org/10.1016/j.jprocont.2017.10.003. Fisher, R. (1935). The design of experiments. Oliver & Boyd, 1935. Ferreau, H.J., Bock, H.G., and Diehl, M. (2008). An online active set strategy to overcome the limitations of explicit MPC. International Journal of Robust and Nonlinear Control, 18(8):816–830, 2008. Heirung, T.A.N., Ydstie, B.E., and Foss, B. (2012). Towards Dual MPC. In Proceedings of the 4th IFAC Nonlinear Model Predictive Control Conference, 45(17): 502–507, 2012. Heirung, T.A.N., Foss, B., and Ydstie, B.E. (2015). MPCbased dual control with online experiment design. Journal of Process Control, 32: 64–76, 2015. Hernandez, B. and Trodden, P. (2015). Persistently exciting tube MPC. American Control Conference (ACC), IEEE, 2016. 517

523

Hovd, M. and Bitmead, R.R. (2005). Interaction between control and state estimation in nonlinear MPC. Modeling, Identification, and Control, 26(3):165–174, 2005. Houska, B., Ferreau, H.J., and Diehl, M. (2011). An autogenerated real-time iteration algorithm for nonlinear MPC in the microsecond range. Automatica, 47(10): 2279-2285, 2011. Houska, B., Frasch, J., and Diehl, M. (2016). An Augmented Lagrangian Based Algorithm for Distributed Non-Convex Optimization. SIAM Journal on Optimization, 26(2):1101–1127, 2016. Houska, B., Telen, D., Logist, F., and Van Impe, J. (2017). Self-reflective model predictive control. SIAM Journal on Control and Optimization, 55(5): 2959–2980, 2017. Kawohl, M., Heine, T., and King, R. (2007). A new approach for robust model predictive control of biological production processes. Chemical Engineering Science, 62(18–20): 5212-5215, 2007. Ljung, L. (1999). System Identification: Theory for the User. Prentice Hall, 1999. Marafioti, G., Bitmead, R.R., and Hovd, M. (2014). Persistently exciting model predictive control. International Journal of Adaptive Control and Signal Processing, 28(6): 536–552, 2014. Mayne, D.Q., Rawlings, J.B., Rao, C.V., and Scokaert, P.Q. (2000). Constrained model predictive control: Stability and optimality. Automatica, 36(6): 789–814, 2000. Mesbah, A., Bombois, X., Forgione, M., Hjalmarsson, H., and Van den Hof, P.M.J. (2015). Least costly closedloop performance diagnosis and plant re-identification. International Journal of Control, 88(11): 2264–2276, 2015. Pukelsheim, F. (1993). Optimal design of Experiments. John Wiley & Sons, Inc., New York, 1993. Rao, C.V., Rawlings, J.B., and Lee, J.H. (2001). Constrained linear state estimation - a moving horizon approach Automatica, 37:1619–1628, 2001. Rawlings, J.B. and Mayne, D.Q. (2009). Model Predictive Control: Theory and Design. Madison, WI: Nob Hill Publishing, 2009. Stengel, R. (1994). Optimal Control and Estimation. Dover Publications, New York, 1994. Telen, D., Houska, B., Logist, F., Van Derlinden, E., Diehl, M., and Van Impe, J. (2013). Optimal Experiment Design under Process Noise. Journal of Process Control, 23: 613–629, 2013. Telen, D., Houska, B., Vallerio, M., Logist, F., and Van Impe, J. (2016). A study of integrated experiment design for NMPC applied to the Droop model. Chemical Engineering Science, 160: 370–383, 2016. Yan, J. and Bitmead, R.R. (2005). Incorporating state estimation into model predictive control and its application to network traffic control. Automatica, 41(4): 595– 604, 2005. Zacekova, E., Privara, S., and Pcolka, M. (2013). Persistent excitation condition within the dual control framework. Journal of Process Control, 23(9): 1270–1280, 2013. Zhu, Y. (2009). System identification for process control: recent experience and outlook. International Journal of Modelling, Identification and Control, 6(2): 89–103, 2009.