International Journal of Automation and Computing  2018, Vol. 15 Issue (5): 616-624 PDF
A Nash Game Approach to Mixed H2/H Model Predictive Control: Part 3-Output Feedback Case
Department of Chemical Engineering, Indian Institute of Technology Guwahati, Assam 781039, India
Abstract: In this paper, the state-feedback Nash game based mixed H2/H design[1, 2] has been extended for output feedback case. The algorithm is applied to control bioreactor system with a Laguerre-Wavelet Network (LWN)[3, 4] model of the bioreactor. This is achieved by using the LWN model as a deviation model and by successively linearising the deviation model along the state trajectory. For reducing the approximation error and to improve the controller performance, symbolic derivation algorithm, viz., automatic differentiation is employed. A cautionary note is also given on the fragility of the output feedback mixed H2/H model predictive controller[4, 5] due to its sensitivity to its own parametric changes.
Key words: Robust model predictive control     mixed H2/H control     Nash game     output feedback model predictive control (MPC)     automatic differentiation     fragility of controller
1 Introduction

Model predictive control (MPC) algorithms require full knowledge of the system states while using state feedback framework. However, it is not possible to control a process directly using state feedback technique when the states of the process are not available for measurement. In such cases, the controller has to be designed with output feedback framework, i.e., to control the process using the output measurements of the process. Output feedback control is accomplished using either direct output-to-input mapping or using the state feedback control by estimating the states from the process output measurement.

Linear systems, not affected by any disturbance, can be controlled by linear control techniques that employ an observer/estimator. In such cases separation principle holds good and stability can be ensured. However, stability of the closed-loop cannot be ensured, in general, by simply combining a stable estimator with a stable state feedback controller, when the system or its controller is nonlinear (as is the case with MPC for constrained systems) and the state and output disturbances are present[6]. To overcome this drawback some special care has to be taken in the robust controller design to deal with the estimation error along with modeling error.

In this paper, the state-feedback Nash game based mixed MPC (NGM-MPC)[1, 2] design has been extended for output feedback case. However, the extension is not straightforward, as it is an observer based controller.

2 Problem formulation

Consider the dynamic system

 $x_{k+1} =Ax_{k}+B_{0}w_{k}^{0}+B_{1}w_{k}+B_{2}u_{k}$ (1)
 $\tilde{z}_{k} =C_{0}x_{k}+D_{01}w_{k}+D_{02}u_{k}$ (2)
 $z_{k} =C_{1}x_{k}+D_{10}w_{k}^{0}+D_{11}w_{k}+D_{12}u_{k}$ (3)
 $y_{k} =C_{2}x_{k}+D_{20}w_{k}^{0}+D_{21}w_{k}$ (4)

where $\tilde{z}_{k}$ and $z_{k}$ are the controlled output corresponding to ${H}_{2}$ and ${H}_{\infty }$ performance measures, $w_{k}^{0}$ is the stochastic/random process noise (e.g., due to high-gain sensors), $w_{k}$ is the unknown-but-bounded disturbance (e.g., the state estimation error), $y_{k}$ is process measured output - satisfying the following assumptions:

Assumption 1. $(A, B_{2})$ is stabilisable and $(C_{2}, A)$ is detectable.

Assumption 2. $w_{k}\in {P}\subset \ell _{2}$ is a bounded power signal1 and $w_{k}^{0}$ is a normalized Gaussian white noise with zero mean and unit variance. Moreover $w_{k}^{0}\perp w_{j}, \forall k\geq j$.

1 The set of sequences ${{P}}\triangleq x:\| x\|_{{{P}} } < \infty$ is called the set of signals with bounded power. It should be noted that $\|\cdot\|_{{P}}$ is not necessarily a norm, since all the signals in $\ell_2\subset {P}$ have zero power. We assume $w\in{ {P}}$.

Assumption 3. The initial condition $x_{0}\sim {N}(\bar{x}_{0}, R_{0})$2 is independent of $w_{k}^{0}, \forall k\geq 0$.

2 $\bar{x}_0$ represents the mathematical mean value of $x_0$ and $R_0$ represents its variance of the normal distribution.

In ${H}_{\infty }$ control problem, the ${H}_{\infty }$ norm $\Vert T_{z_{k}, w_{k}}\Vert _{\infty }^{2}$ can also be given in terms of the semi-norm $\Vert \cdot \Vert _{{P}}$ for a disturbance bound $\gamma$ as

 \begin{align} \Vert T_{z_{k}, w_{k}}\Vert _{\infty }^{2}=\sup\limits_{w_{k}\in {P}}\frac{ \Vert z_{k}\Vert _{{P}}^{2}}{\Vert w_{k}\Vert _{{P}}^{2}} <\gamma ^{2} \end{align} (5)

whose equivalent form is

 \begin{align} \sup\limits_{w_{k}\in {P}}\{\Vert z_{k}\Vert _{{P}}^{2}-\gamma ^{2}\Vert w_{k}\Vert _{{P}}^{2}\}. \end{align} (6)

Being a multi-objective control problem, where the two performance functions, namely, ${H}_{2}$ and ${H}_{\infty }$, are given, respectively, as follows:

 $J_{2}(w_{k}, u_{k})= \arg \min\limits_{u\in {U}}\Vert T_{\tilde{z} _{k}, w_{k}^{0}}\Vert _{2}^{2}$ (7)
 $J_{\infty }(w_{k}, u_{k})= \sup\limits_{w_{k}\in {P}}\{\Vert z_{k}\Vert _{ {P}}^{2}-\gamma ^{2}\Vert w_{k}\Vert _{{P}}^{2}\}$ (8)

where $T_{\tilde{z}_{k}, w_{k}^{0}}$ represents the transfer function from $w_{k}^{0}$ to $\tilde{z}_{k}$.

As the state feedback problem is solvable, the output feedback problem is also solvable once the states are estimated by an observer. Then the resulting controller will be given as

 \begin{align} u^*_k=K_{k}\hat{x}_k \end{align} (9)

where $\hat{x}_k$ is the estimated state of the observer.

2.1 Solving the ${\pmb H_2}$ problem

Let $w_{k}^{\ast }=H_{k}x_{k}$ be the worst case disturbance for the full information problem. On substituting this in the original system equation, one can safely eliminate the ${H}_{\infty }$ part from the mixed ${H}_{2}$/${H}_{\infty }$ problem, i.e., (3) will not be considered in this part of solving the ${H} _{2}$ problem.

 $x_{k+1} = (A+B_{1}H_{k})x_{k}+B_{0}w_{k}^{0}+B_{2}u_{k}$ (10)
 $\tilde{z}_{k} =(C_{0}+D_{01}H_{k})x_{k}+D_{02}u_{k}$ (11)
 $y_{k} = (C_{2}+D_{21}H_{k})x_{k}+D_{20}w_{k}^{0}.$ (12)

Now the problem to be solved is

 \begin{align} u_{k}^{\ast }=\arg \min\limits_{u}\parallel T_{\tilde{z}_{k}, w_{k}^{0}}\parallel _{2}^{2}. \end{align} (13)

Generally, for a discrete time linear system given by

 $x_{k+1} = \bar{A}x_{k}+\bar{B}u_{k}+\bar{E}w_{k}^{0}$ (14)
 $y_{k} = \bar{C}_{1}x_{k}+\bar{D}_{1}w_{k}^{0}$ (15)
 $z_{k} = \bar{C}_{2}x_{k}+\bar{D}_{2}u_{k}$ (16)

the solution of discrete-time ${H}_{2}$ optimal control problem is given by the following two discrete-time algebraic Riccati equations (DAREs)[7]

 \begin{align} P =&\bar{A}^{\rm T}P\bar{A}+\bar{C}_{2}^{\rm T}\bar{C}_{2}-(\bar{C}_{2}^{\rm T}\bar{D} _{2}+\bar{A}^{\rm T}P\bar{B})\times \notag \\ &(\bar{D}_{2}^{\rm T}\bar{D}_{2}+\bar{B}^{\rm T}P\bar{B})^{-1}(\bar{D}_{2}^{\rm T}\bar{C }_{2}+\bar{B}^{\rm T}P\bar{A}) \end{align} (17)
 \begin{align} Q =&\bar{\bar{A}}Q\bar{\bar{A}}^{\rm T}+\bar{E}\bar{E}^{\rm T}-(\bar{\bar{A}}Q\bar{C }_{1}^{\rm T}+\bar{E}\bar{D}_{1}^{\rm T})\times \notag \\ &(\bar{D}_{1}\bar{D}_{1}^{\rm T}+\bar{C}_{1}Q\bar{C}_{1}^{\rm T})^{-1}(\bar{D}_{1} \bar{E}^{\rm T}+\bar{C}_{1}Q\bar{\bar{A}}^{\rm T}) \end{align} (18)

whose brief derivation is given in Appendix A. The above pair of DAREs could also be rewritten using matrix inversion identity3 as

3 $I_{n\times n}- \tilde{B}(\tilde{C}\tilde{B}+\tilde{D})^{-1}\tilde{C}=(I_{n\times n}+\tilde{B }\tilde{D}^{-1}\tilde{C})^{-1}$, where $\tilde{B}\in {\bf R}^{n\times 1}$, $\tilde{C}\in {\bf R}^{1\times n}$, $\tilde{D}\in {\bf R}^{1\times 1}$.

 \begin{align} P =&\bar{A}^{\rm T}P\left\{ I+\bar{B}(\bar{D}_{2}^{\rm T}\bar{D}_{2})^{-1}\bar{B} ^{\rm T}P\right\} ^{-1}\bar{A} + \nonumber \\ &\bar{C}_{2}^{\rm T}\left\{ I-\bar{D}_{2}^{\rm T}(\bar{D}_{2}\bar{D}_{2}^{\rm T})^{-1} \bar{D}_{2}^{\rm T}\right\} \bar{C}_{2} \end{align} (19)
 \begin{align} Q =&\bar{\bar{A}}Q\left\{ I-\bar{C}_{1}^{\rm T}(\bar{D}_{1}\bar{D}_{1}^{\rm T})^{-1} \bar{C}_{1}Q\right\} ^{-1}\bar{\bar{A}} + \nonumber \\ &\bar{E}\left\{ I-\bar{D}_{1}^{\rm T}(\bar{D}_{1}\bar{D}_{1}^{\rm T})^{-1}\bar{D} _{1}\right\} \bar{E}^{\rm T}. \end{align} (20)

Using the analogy of the above general ${H}_{2}$ control problem in (14)-(16) with the given system information in (10)-(12) and by assigning

 $\bar{A} =A_{C}+(B_{1}-B_{2}R_{02}^{-1}D_{02}D_{01})H_{k}$ (21)
 $\bar{\bar{A}} =A_{O}+(B_{1}-B_{0}D_{20}^{\rm T}R_{20}^{-1}D_{21})H_{k}$ (22)
 $\bar{B} =B_{2}$ (23)
 $\bar{C} =C_{0}+D_{21}H_{k}$ (24)
 $\bar{E} =B_{0}$ (25)
 $\tilde{R}_{02}= I-\bar{D}_{2}^{\rm T}(\bar{D}_{2}\bar{D}_{2}^{\rm T})^{-1}\bar{D}_{1}$ (26)
 $\tilde{R}_{20}= I-\bar{D}_{1}^{\rm T}(\bar{D}_{1}\bar{D}_{1}^{\rm T})^{-1}\bar{D}_{1}.$ (27)

we can get the solution for the ${H}_{2}$ part of the mixed ${H}_{2}$/${H}_{\infty }$ problem by solving the following two algebraic Riccati equations

 \begin{align} Ric1 :&0=\left\{ \begin{array}{l} \left[ A_{C}+\left( B_{1}-B_{2}R_{02}^{-1}D_{02}D_{01}\right) H_{k}\right] ^{\rm T}\times \\ P_{k}^{1}\left( I+B_{2}R_{02}^{-1}B_{2}^{\rm T}P_{k}^{1}\right) ^{-1}\times \\ \left[ A_{C}+\left( B_{1}-B_{2}R_{02}^{-1}D_{02}D_{01}\right) H_{k}\right] -\\ P_{k}^{1}+\left( C_{0}+D_{21}H_{k}\right) ^{\rm T}\tilde{R}_{02}\times \\ \left( C_{0}+D_{21}H_{k}\right) \end{array} \right. \label{Ric1Out} \end{align} (28)
 \begin{align} Ric2 :&0=\left\{ \begin{array}{l} \left[ A_{0}+\left( B_{1}-B_{0}D_{20}^{\rm T}R_{20}^{-1}D_{21}\right) H_{k} \right] \times \\ P_{k}^{2}\left( I-C_{2}^{\rm T}R_{20}^{-1}C_{2}P_{k}^{2}\right) ^{-1}\times \\ \left[ A_{0}+\left( B_{1}-B_{0}D_{20}^{\rm T}R_{20}^{-1}D_{21}\right) H_{k} \right] ^{\rm T} -\\ P_{k}^{2}+B_{0}\tilde{R}_{20}B_{0}^{\rm T} \end{array} \right. \label{Ric2Out} \end{align} (29)

where

 $A_{C} =A-B_{2}R_{02}^{-1}D_{02}^{\rm T}C_{0}$ (30)
 $A_{O} =A-B_{0}D_{02}^{\rm T}R_{20}^{\rm T}C_{2}$ (31)
 $R_{02} =D_{02}^{\rm T}D_{02}$ (32)
 $\tilde{R}_{02} =I-D_{02}R_{02}^{-1}D_{02}^{\rm T}$ (33)
 $R_{20} =D_{20}D_{20}^{\rm T}$ (34)
 $\tilde{R}_{20} =I-D_{20}R_{20}^{-1}D_{20}$ (35)

and the state feedback controller and observer gains are

 \begin{align} K_{k} =&-\left[ D_{02}^{\rm T}D_{02}+B_{2}^{\rm T}P_{k}^{1}B_{2}\right] ^{-1}\times \nonumber \\ & \left[ \begin{array}{l} B_{2}^{\rm T}P_{k}^{1}(A+B_{1}H_{k})+ \\ D_{02}^{\rm T}(C_{0}+D_{01}H_{k}) \end{array} \right] \end{align} (36)
 \begin{align} L_{k} =&-[(A+B_{1}H_{k})P_{k}^{2}(C_{2}+D_{21}H_{k})^{\rm T}+B_{0}D_{20}^{\rm T}] \times\nonumber \\ & \left[ \begin{array}{l} D_{20}D_{20}^{\rm T}+(C_{2}+D_{21}H_{k})P_{k}^{2}\times \\ (C_{2}+D_{21}H_{k})^{\rm T} \end{array} \right] ^{-1} \end{align} (37)

where $P_{k}^{1}$ and $P_{k}^{2}$ are the positive semi-definite matrices giving the stabilising solution to the ${H}_{2}$ problem.

This part of the control problem finds out the optimal controller such that $y\mapsto u$, which has the realization given by

 \begin{align} \tilde{x}_{k+1} =&\left( \begin{array}{c} A+B_{1}H_{k}+B_{2}K_{k}+ \\ L_{k}C_{2}+L_{k}D_{21}H_{k} \end{array} \right) \tilde{x}_{k}-L_{k}y_{k} \end{align} (38)
 \begin{align} u_{k}^{\ast } =&K_{k}\tilde{x}_{k}. \end{align} (39)

The above two equations represent the output-feedback controller, the former is an observer and the latter is a state feedback controller.

2.2 Solving the ${\pmb H_{\infty}}$ problem

Now, with the optimal control law $u_{k}^{\ast }=K_{k}\hat{x}_{k}$, the performance function $J_{\infty }(w, u)$ for the system given in (1) and (3), is given as

 $x_{k+1} =(A+B_{2}K_{k})x_{k}+B_{0}w_{k}^{0}+B_{1}w_{k}$ (40)
 $z_{k} =(C_{1}+D_{12}K_{k})x_{k}+D_{10}w_{k}^{0}+D_{11}w_{k}.$ (41)

Using the Stochastic Bounded Real Lemma[8] for a general discrete-time system

 $x_{k+1} =\bar{A}x_{k}+\bar{B}_{0}w_{k}^{0}+\bar{B}_{1}w_{k}$ (42)
 $z_{k} =\bar{C}_{1}x_{k}+\bar{D}_{10}w_{k}^{0}+\bar{D}_{11}w_{k}$ (43)

the corresponding general discrete-time ${H}_{\infty }$ algebraic Riccati equation is given as

 \begin{align} 0 =&\bar{A}^{\rm T}X\bar{A}-X+(\bar{B}_{1}^{\rm T}X\bar{A}+\bar{D}_{11}^{\rm T}\bar{C} _{1})^{\rm T} \times \notag \\ & (\gamma ^{2}I-\bar{D}_{11}^{\rm T}\bar{D}_{11}-\bar{B}_{1}^{\rm T}X\bar{B} _{1})^{-1}\times \notag \\ & (\bar{B}_{1}^{\rm T}X\bar{A}+\bar{D}_{11}^{\rm T}\bar{C}_{1})+\bar{C} _{1}^{\rm T}\bar{C}_{1} \label{Ric_22} \end{align} (44)

where $X\succ 0$. For some $\gamma > 0$, the worst case bounded disturbance signal is given by

 \begin{align} w_{k}^{\ast }=(\gamma ^{2}I-\bar{D}_{11}^{\rm T}\bar{D}_{11}-\bar{B}_{1}^{\rm T}X \bar{B}_{1})^{-1}(\bar{B}_{1}^{\rm T}X\bar{A}+\bar{D}_{11}^{\rm T}\bar{C}_{1})x_{k}. \label{uCons} \end{align} (45)

In the light of the analogy of structure of the above given general discrete-time system with the system given by (40)-(41), one can get

 $\bar{A} =A+B_{2}K_{k}$ (46)
 $\bar{B}_{0} =B_{0}$ (47)
 $\bar{B}_{1} =B_{1}$ (48)
 $\bar{C}_{1} =C_{1}+D_{12}K_k$ (49)
 $\bar{D}_{10} =D_{10}$ (50)
 $\bar{D}_{11} =D_{11}$ (51)

and the corresponding discrete-time ${H}_{\infty }$ algebraic Riccati equation as

 \begin{align} Ric3:\, 0=\left\{ \begin{array}{l} \left( A+B_{2}K_{k}\right) ^{\rm T}P_{k}^{3}\left( A+B_{2}K_{k}\right) -P_{k}^{3}+ \\ \left[ B_{1}^{\rm T}P_{k}^{3}\left( A+B_{2}K_{k}\right) +D_{11}^{\rm T}\left( C_{1}+D_{12}K_{k}\right) \right] ^{\rm T} \\ \left( \gamma ^{2}I-D_{11}^{\rm T}D_{11}-B_{1}^{\rm T}P_{k}^{3}B_{1}\right) ^{-1} \\ \left[ B_{1}^{\rm T}P_{k}^{3}\left( A+B_{2}K_{k}\right) +D_{11}^{\rm T}\left( C_{1}+D_{12}K_{k}\right) \right]+ \\ \left( C_{1}+D_{12}K_{k}\right) ^{\rm T}\left( C_{1}+D_{12}K_{k}\right) \end{array} \right. \label{Ric3Out} \end{align} (52)

where $P_{k}^{3}\succ 0$ is the stabilizing solution of the algebraic Riccati equation, with the worst-case disturbance gain

 \begin{align} H_{k}=(\gamma ^{2}I-D_{11}^{\rm T}D_{11}-B_{1}^{\rm T}P_{k}^{3}B_{1})^{-1}(B_{1}^{\rm T}P_{k}^{3}A+D_{11}^{\rm T}C_{1}) \end{align} (53)

and

 \begin{align} w_{k}^{\ast }=H_{k}\hat{x}_{k}. \end{align} (54)

Thus (28), (29) and (53) are the three coupled algebraic Riccati equations (3-cAREs) of the optimal output feedback control problem. The coupling is attributed due to the presence of the ${H}_{\infty }$ solution term, $H_{k}$, in (28) and (29); similarly, the solution of the ${H}_{2}$ part of the problem, $K_{k}$, in the (52).

The saddle point solution, viz. Nash equilibria, for the output feedback case is similar to that of the state feedback case, except that the state information is been estimated from the observer:

 $\frac{1}{2}\hat{x}_{k}^{\rm T}P_{k}^{1}\hat{x}_{k} \geq \beta _{k}$ (55)
 $\frac{1}{2}\hat{x}_{k}^{\rm T}P_{k}^{3}\hat{x}_{k} \leq \beta _{k}.$ (56)

With this bound, updating the control law using the current state of the system, the optimal values of the stabilizing matrices viz. $P_{k}^{1}$, $P_{k}^{2}$ and $P_{k}^{3}$ are found out.

3 Solving the 3-cAREs

This is again similar to the case of state feedback control problem. However, it involves 3-cAREs in the case of output feedback optimal control problem.

Theorem 1. For a system defined by $\Sigma :(A, B_{0}, B_{1}, B_{2})$ with initial condition $x_{0}$, the two-player Nash game based mixed $H_{2}$/$H_{\infty }$ MPC yields an (sub-)optimal, stable solution, for the three cross-coupled algebraic Riccati equations (3-cAREs)

 $Ric1+\epsilon _{k}^{1} =0$ (57)
 $Ric2+\epsilon _{k}^{2} =0$ (58)
 $Ric3+\epsilon _{k}^{3} =0$ (59)

and along with the saddle point optimal condition bounds

 $\hat{x}_{k}^{\rm T}P_{k}^{1}\hat{x}_{k} \geq \beta _{k}$ (60)
 $\hat{x}_{k}^{\rm T}P_{k}^{3}\hat{x}_{k} \leq \beta _{k}$ (61)

if it exists, such that ($P_{k}^{1^{\ast }}, P_{k}^{2^{\ast }}, P_{k}^{3^{\ast }})\succ 0$, for some open-hyper balls of radius $\left\{ \epsilon _{k}^{1^{\ast }}, \epsilon _{k}^{2^{\ast }}, \epsilon _{k}^{3^{\ast }}\right\} \in {\bf R}^{n\times n}$ around and containing the origin, such that the closed loop system is stable.

Proof. Let $P_{k}^{1^{\ast }}, P_{k}^{2^{\ast }}$ and $P_{k}^{3^{\ast }}$ be the (sub-)optimal, stable and locally unique solution set of the 3-cAREs (28, 29 and 52), satisfying the conditions (60 and 61) and $u_{k}\in {U}$ (i.e., satisfying 45), with $x_{0}\in {X}_{0}$, for some open hyper balls of radius $\left\{ \epsilon _{k}^{1^{\ast }}, \epsilon _{k}^{2^{\ast }}, \epsilon _{k}^{3^{\ast }}\right\} \in {\bf R}^{n\times n}$, respectively, in the linear space. Now assume that there exists another solution set, $\tilde{P}_{k}^{1\ast }$, $\tilde{P}_{k}^{2\ast }$ and $\tilde{P}_{k}^{3\ast }$, which again satisfies (28), (29) and (52) for some other hyper balls of radius $\{\tilde{\epsilon}_{k}^{1^{\ast }}, \tilde{\epsilon} _{k}^{2^{\ast }}, \tilde{\epsilon}_{k}^{3^{\ast }}\}\in {\bf R}^{n\times n}$, respectively, such that $\epsilon _{k}^{j^{\ast }}=\tilde{\epsilon} _{k}^{j^{\ast }}+\delta _{j}\in {\bf R}^{n\times n}$ and $\delta _{j}\ll \tilde{\epsilon}_{k}^{j^{\ast}}, \; j=1, 2, 3$.

Likewise, there may exist a different solution set for a different choice of $\delta _{j}$ for every other choice of solution set. Therefore, as $|\delta _{j}|\rightarrow 0$, in the linear space the open hyper balls converges to the (sub-)optimal solution

 \begin{align} \lim\limits_{|\delta _{j}|\rightarrow 0}\epsilon _{k}^{j^{\ast }}=\tilde{\epsilon} _{k}^{j^{\ast }}+\delta _{j}=\tilde{\epsilon}_{k}^{j^{\ast }}, \;\;j=1, 2, 3 \end{align} (62)

which in turn gives $P_{k}^{1^{\ast }}=\tilde{P}_{k}^{1^{\ast }}$, $P_{k}^{2^{\ast }}=\tilde{P}_{k}^{2^{\ast }}$ and $P_{k}^{3^{\ast }}=\tilde{P} _{k}^{3^{\ast }}$, thus giving an (sub-)optimal, stable and locally unique solution set to the 3-cAREs.

One of the main objectives of this paper is to extend the mixed ${H} _2$/${H}_{\infty}$ MPC for the output feedback case to use it with the Laguerre-Wavelet network (LWN) model. For doing so, as the static nonlinear wavelet network portion cannot be inverted, it has to be linearised online, using some efficient numerical techniques such as automatic differentiation (AD). The use of tools such as AD for linearisation not only reduces the computational demand of numerical differentiation but also gives error-free differentiation values, usually up to the machine precision level.

4 Successive linearisation model

Model predictive control relies explicitly on the process model to make an optimal decision on the control input/law to be implemented at every time instant. Hence the correctness of the process model at the prevailing operating condition, is a major factor for efficient performance of MPC. The Wiener type Laguerre-wavelet network model[3] is a data-driven nonlinear model. Usually the input-output data collected for identifying nominal stable systems are obtained by exciting the process about its nominal equilibrium point. Hence the efficient control of the process by linearising the model at the nominal operating point and controlling the process using linear controller is possible only in the close vicinity of the nominal operating point of the process. However, in situations where the process operating point is shifted far away from the previous value, the above mentioned model (linearised at the nominal operating point)-based control strategy can fail miserably.

Successive linearisation approach[9] could be adopted to circumvent the above shortcoming. The successive linearisation of the static nonlinear gain by finding the Jacobian linearisation of the nonlinearity around the nominal state trajectory is made to obtain a local linear equivalent model of the state-to-output mapping. So linearising the wavelet network about the nominal value of the current Laguerre states gives a local linear state-to-output gain. The necessary condition to ensure a smooth linearized representation is that the static nonlinear gain $\Psi$[3, 10] should be Lipschitz continuous. Thus instead of a priori offline piecewise linearisation over the entire range of the nonlinear gain and storing the models to represent the system as a piecewise affine system, the online successive linearisation approach demands lesser memory requirement. Although the associated computational cost and the approximation error may be forbidding with earlier numerical techniques to find the Jacobian, the numerical derivative techniques such as automatic differentiation[11-13] too, could be used to overcome these issues. So automatic differentiation method is used for the successive linearisation using the Matlab automatic differentiation (MAD) module of TOMLAB version 6.1 in the present work.

4.1 Automatic differentiation

From the standard reference for the subject[11] it is stated that:

Algorithmic or automatic differentiation (AD) is concerned with the accurate and efficient evaluation of derivatives for functions defined by computer programs.

AD is a numerical differentiation technique that uses chain-rule of calculus (differentiation) for the floating point evaluation of a function and/or its derivatives[13]. The principal difference between AD and the usual numerical finite difference methods is that in AD there is no discretisation or cancellation errors; moreover, the results (derivative values) are accurate to the machine precision or within the round-off used. Symbolic differentiation methods makes a function derivative into a single long expression at the point where the function$'$s derivative has to be evaluated. This drastically increases the memory requirement for evaluating the derivatives. However, on the other hand, AD uses the common control structures such as loops, branches, sub-functions, etc., unlike symbolic differentiation methods. Also, AD works on floating point numbers making the results more accurate. Thus AD is more efficient than symbolic differentiation.

Usually AD is implemented in either of the two ways: 1) operator overloading, 2) source transformation. In operator overloading the existing classes or types are redefined, initialised with appropriate values for the function and its derivative, invoke the function, and find the values of the derivatives. Source transformation is a method that requires suitable sophisticated compiler-like software to read the computer program containing the function and determine which statements for the function$'$s program requires differentiation. Accordingly, a new version of the given/original program is generated augmenting with the statements for which derivatives need to be calculated.

There are fundamentally two approaches for computing the derivatives using AD.

1) Forward mode: The function is evaluated for its derivative in the forward direction starting from the input (or independent) variable. In this forward mode of calculation, the sensitivities for finding the derivatives of the intermediate stages are propagated forward until the final result.

2) Reverse mode: This works in two stages. In the first stage the original code is run and is augmented with the statements to store data. Then the derivative is evaluated starting from the output (dependent variable) towards the input using all the sensitivities (also called adjoints) in the reverse direction in its second run.

In the present work, TOMLAB-MAD has been used for the purpose of linearisation, so as to reduce the error, incurred due to the evaluation of derivative values computationally, by methods such as finite difference. The use of AD for finding the gradient or Jacobian using MAD software package is efficient and fast enough to implement it online with a regular PC.

5 Numerical example

To demonstrate the performance of mixed ${H}_{2}$/${H} _{\infty }$ output feedback MPC, a simple stable discrete time LTI system is taken. Choice of one such example could be justified with the claim that the developed controller will be used with the proposed Laguerre wavelet network model, where the linear dynamic part of the of the model, i.e., Laguerre filter will always have stable poles. A second order stable system, shown below, with initial condition $x_{0}=[1\; 0]^{\rm T}$ is taken for this purpose. The minimal workable $\gamma$ value for the present example using the proposed robust MPC algorithm is 5.

 \begin{align*} A =&\left[ \begin{array}{cc} 0.80&-0.09 \\ 1.00&0.00 \end{array} \right] \end{align*}
 \begin{align*} B_{0} =&\left[ \begin{array}{c} 0.1 \\ 1.0 \end{array} \right] , B_{1}=\left[ \begin{array}{c} 0.1 \\ 0.0 \end{array} \right], B_{2}=\left[ \begin{array}{c} 0.0 \\ 1.0 \end{array} \right]\\ C_{0} =&\left[ \begin{array}{ll} 1&1 \end{array} \right] , C_{1}=\left[ \begin{array}{ll} 1&0 \end{array} \right] , C_{2}=\left[ \begin{array}{ll} 1.1&-0.07 \end{array} \right] \\ D_{00} =&1.0, D_{01}=0.1, D_{02}=0.001, \\ D_{10} =&0.01, D_{11}=1.0, D_{12}=0.1, \\ D_{20} =&0.1, D_{21}=0.0. \end{align*}

The system is perturbed with the deterministic disturbance, $w_{k}=0.1$ and the stochastic Gaussian white noise ($w_{k}^{0}$) has mean value of zero and standard deviation as 0.1 for all $k\geq 0$. The closed loop performance is given in Figs. 1 and 2.

 Download: larger image Fig. 1. Input and output profiles of the Nash game based mixed ${H}_{2}$/${H}_{\infty }$ MPC

 Download: larger image Fig. 2. State estimation error profiles of the Nash game based mixed ${H}_{2}$/${H}_{\infty }$ MPC

It can be observed that the effect of the measurement noise component $w_{k}^{0}$ does not affect that controller stability. Moreover, the presence of persistent bounded power signal $w_{k}$ has its effect on the state estimation, which does not allow the estimation error of one of states ($e_{2}$) to decay to zero. However, the presence of such a non-zero bounded power disturbance does not affect the closed-loop stability of the process.

6 Bioreactor control

Control of bioreactor is of much interest due to its highly nonlinear behaviour. A mechanistic model[14] of the bioreactor is given in the Appendix B. The dynamics of bioreactor, which exhibits input-multiplicity nonlinearity, is shown elsewhere[14, 15]. At the low concentration of input feed substrate, the growth rate of the cell-biomass dominates the reaction rate. Thus, increase in feed concentration results in an increase in productivity. However, beyond a certain value of the input feed concentration, the product and substrate inhibition prevails and thereupon productivity shows a negative gradient with respect to the input feed concentration. Thus, there exists an optimum value of the bioreactor productivity. The necessary condition for optimality implies that the steady-state gain (product concentration vs. manipulated feed substrate concentration) changes its sign across the optimum point and the steady state gain is zero at the optimum. Thus, the system exhibits input multiplicity and there exists a singularity at the optimum operating point.

The robust control of the bioreactor is one of the objectives of the paper. For achieving this, the dynamics of bioreactor system is efficiently captured through the nonlinear system identification method using LWN model from the process input-output data[3]. The mechanistic model (Appendix B) of the bioreactor[16] is assumed to be the actual process for the simulation case studies.

6.1 Control of bioreactor using mixed ${\pmb {H_2/H _{\infty}}}$ output feedback MPC

In the present study on the bioreactor process, the closed loop robustness due to mixed ${H}_{2}$/${H}_{\infty }$ output feedback MPC is manifested against system parametric uncertainty. For any data-driven model developed using a suitable system identification technique, e.g., the LWN model for the bioreactor process under consideration, such discrepancies eventually lead to plant-model mismatch. Had there been any adaptive mechanism as adopted by [16], this plant-model mismatch can be circumvented in an efficient way. However, this adaptive procedure is not always feasible in real-time application mainly due to the complexity of the model and/or the computational cost associated with it. In such cases, the burden comes to the controller design, to achieve the desired plant operation. In the present bioreactor example, the nonlinearity of the process is due to state-to-output nonlinear static gain, which is represented by the wavelet network in the Wiener type LWN model.

For some processes the plant-model mismatch may result in deterioration of the performance or may even lead to catastrophic effects, depending on the nature of the process. The state-feedback controller design[3, 10] cannot be deployed for such processes, where the system states are not directly measurable for control, or, for model based control cases, where the states of the model do not reflect the actual process states. For any black-box model, usually, the model states fall under the latter category. In such cases, it is wise to use successively linearized deviation (SLD) model, than simply employing successive linearisation technique. State-deviation (SD) model is one in which the model, whose states represent deviation of the states from a reference value (usually the previous value of the model state) as $\tilde{x}_{k}:=x_{k}-x_{k-1}$, developed using the identification technique. Then w.r.t. the deviation state, at every time instant, the static nonlinear portion of the Wiener model is linearized. This will result in an SLD model. The SLD model is sought because the NGM-MPC control algorithms (both state-feedback and output feedback) consider the origin as the nominal equilibrium point, i.e., $\{0, 0\}\in \text{int }\{{X}\times {U}\}$. The response of the SLD-LWN model for the bioreactor process, excited by a sinusoidal input signal about the nominal input feed concentration of 20 g/L is shown in Fig. 3.

 Download: larger image Fig. 3. SLD-LWN model response for sinusoidal input signal

It could be observed that there is not much distortion in the performance of the LWN model due to successively linearised deviation approach. Hence the approximation error of the nominal model due to the linearisation approach is minimal and thereby the controller$'$s capacity is maximally exploited to handle plant-model mismatch, due to approximation or worst-case disturbance. The plant-model mismatch can occur due to the parametric change in the process. For instance, the maximum specific growth rate of the microbes ($\mu_{m}$) is taken as 0.48 in the present nominal case. The nominal LWN model is developed using the input-output data of the process only for this value of $\mu_{m}$. However, there could always exist the possibilities of change in the maximum specific growth of the microbes ($\mu_{m}$) due to external reasons such as change in ambient temperature, etc. So the LWN model developed using the nominal value cannot represent the actual process in those circumstances. This will cause a shift in the nominal operating point of the process and hence introduces plant-model mismatch.

The closed-loop performance of the output feedback NGM-MPC controller for step change in $\mu _{m}$ from its nominal value of $0.48$ to $0.44$ and $0.46$ are given in Figs. 4 and 5, respectively. The value of $\gamma$ taken in both the cases for the controller design is 1.55. Change in the process parameter is also considered as disturbance rejection control problem and the process is controlled using the output-feedback NGM-MPC given in the present work.

 Download: larger image Fig. 4. Parametric uncertainty: Output and input profiles of the bioreactor using Nash game based mixed ${H}_{2}$/${H}_{\infty }$ MPC for a step change $\left(\mu _{m}=0.44\right)$

 Download: larger image Fig. 5. Parametric uncertainty: Output and input profiles of the bioreactor using Nash game based mixed ${H}_{2}$/${H}_{\infty }$ MPC for a step change $\left(\mu _{m}=0.46\right)$

6.2 Comparison of output feedback MPC Schemes

The efficacy of the proposed mixed ${H}_{2}$/${H}_{\infty }$ output-feedback MPC (Controller A) is compared with another output feedback MPC design technique (Controller B[17]), for a class of nonlinear systems represented similarly by Wiener type nonlinear model. Although the model considered by [17] is Hammerstein-wiener type, the Hammerstein portion at the input side of the model is removed by means of a suitable inverse function. The static nonlinearity at the output side of the model is interestingly represented by suitable linear inclusion as polytopic descriptions for controlling a class of constrained nonlinear systems. Moreover, the referred paper[17] deals with an infinite horizon min-max output feedback MPC problem to guarantee closed-loop stability. These are some common design similarities of reference[17] with that of the output feedback controller developed in this paper, that makes comparison of their performances reasonable. In this paper, comparison of the closed-loop controller performance for controlling an uncertain nonlinear process viz., continuous bioreactor, with the proposed controller design and that by [17] is made. The comparison is made for a parametric change of $\mu _{m}$ from its nominal value of 0.48 to 0.46 from the 10th hour. Surprisingly, there is hardly any control action by the Controller B given in [17]. However, the proposed Controller A stabilizes the process and bring the process variable to its nominal value despite the parametric uncertainty. Although the failure of Controller B is surprising, the reason for this failure could be attributed to the fact that the linearising technique used by the authors of Controller B are never mentioned explicitly and it is perhaps the fragility of the linearisation that fails the controller. Nevertheless, in this paper, the linearising technique used for controller A & B are the same viz., Automatic differentiation based successive linearisation technique. Thus the effect of the comparison is fairly justified.

7 Conclusions

In this paper, an observer based infinite-horizon output feedback mixed ${H}_{2}$/${H}_{\infty }$ MPC is designed, for those cases where the actual process states are not directly measurable. The control problem results in solving three coupled algebraic Riccati equations. Moreover, as the output feedback MPC is designed for linear dynamic processes, to use it with a class of nonlinear processes represented by Wiener type models[3], successive linearized deviation (SLD) model is used. For reducing the error due to numerical method of linearisation, algorithmic/automatic differentiation technique is advocated. The LWN model developed by [3] is used in the output feedback mixed ${H}_{2}$/${H}_{\infty }$ MPC design for controlling the benchmark process (bioreactor).

Appendix A Discrete-time algebraic Riccati equations of ${\pmb {H_2}}$ optimal control

For a discrete time linear system given by

 \begin{align} x_{k+1} =&\bar{A}x_{k}+\bar{B}u_{k}+\bar{E}w_{k}^0 \end{align} (A1)
 \begin{align} y_{k} =&\bar{C}_{1}x_{k}+\bar{D}_{1}w_{k}^0 \end{align} (A2)
 \begin{align} z_{k} =&\bar{C}_{2}x_{k}+\bar{D}_{2}u_{k} \end{align} (A3)

let the energy (Lyapunov) function associated with it is defined as $V_{k}(x_{k}):=x_{k}^{\rm T}P_{k}x_{k}$.

As the horizon tends towards infinity in an infinite-horizon control problem, it is known that

 \begin{align} \Delta V_k(x_k)=V_k(x_{k+1})-V_k(x_k)=0. \end{align} (A4)

Moreover, the controlled output $z_k$ is assumed to reach their nominal values, so that

 \begin{align} \label{zTz} z_k^{\rm T}z_k=0. \end{align} (A5)

Using the definition of the Lyapunov function and the dynamic state equation, $x_{k+1}=Ax_{k}+Bu_{k}$, for part of the controller design in ${H}_{2}$ optimal control problem, we get

 \begin{align} \Delta V_{k}(x_{k}) =&(\bar{A}x_{k}+\bar{B}u_{k})^{\rm T}P_{k}(\bar{A}x_{k}+ \bar{B}u_{k})-x_{k}^{\rm T}P_{k}x_{k} = \notag \\ &x_{k}^{\rm T}\bar{A}^{\rm T}P\bar{A}x_{k}+x_{k}^{\rm T}\bar{A}^{\rm T}P_{k}\bar{B} u_{k}+u_{k}^{\rm T}\bar{B}^{\rm T}P_k\bar{A}x_k+\notag \\ &u_{k}^{\rm T}\bar{B}^{\rm T}P_{k}\bar{B}u_{k}-x_{k}^{\rm T}P_{k}x_{k}.\label{dVex} \end{align} (A6)

On the other hand, (A5) gives

 \begin{align} z_{k}^{\rm T}z_{k}=[x_{k}^{\rm T}\;\;u_{k}^{\rm T}]\left[ \begin{array}{c} \bar{C}_{2}^{\rm T} \\ \bar{D}_{2}^{\rm T} \end{array} \right] [\bar{C}_{2}\;\;\bar{D}_{2}]\left[ \begin{array}{c} x_{k} \\ u_{k} \end{array} \right] =0. \label{zTzex} \end{align} (A7)

Therefore, the summation of (A6) and (A7) also equates to zero i.e., $S:=\Delta V(x_{k})+z_{k}^{\rm T}z_{k}$

 \begin{align} S &=[x_{k}^{\rm T}\;\;u_{k}^{\rm T}]\left( \begin{array}{c} \left[ \begin{array}{cc} \bar{A}^{\rm T}P\bar{A}-C&A^{\rm T}PB \\ \bar{B}^{\rm T}P\bar{A}&\bar{B}^{\rm T}P\bar{B} \end{array} \right] + \\[4mm] \left[ \begin{array}{cc} \bar{C}_{2}^{\rm T}\bar{C}_{2}&\bar{C}_{2}^{\rm T}\bar{D}_{2} \\ \bar{D}_{2}^{\rm T}\bar{C}_{2}&\bar{D}_{2}^{\rm T}\bar{D}_{2} \end{array} \right] \end{array} \right) \left[ \begin{array}{c} x_{k} \\ u_{k} \end{array} \right] = \notag \\ &\quad [x_{k}^{\rm T}\;\;u_{k}^{\rm T}]\times \notag \\ &\left[ \begin{array}{cc} \bar{A}^{\rm T}P\bar{A}-P+\bar{C}_{2}^{\rm T}\bar{C}_{2}&\bar{A}^{\rm T}P\bar{B}+\bar{C }_{2}^{\rm T}\bar{D}_{2} \\ \bar{B}^{\rm T}P\bar{A}+\bar{D}_{2}^{\rm T}\bar{C}_{2}&\bar{B}^{\rm T}P\bar{B}+\bar{D} _{2}^{\rm T}\bar{D}_{2} \end{array} \right]\times \notag \\ & \left[ \begin{array}{c} x_{k} \\ u_{k} \end{array} \right]=0. \end{align} (A8)

From the above, we get the following algebraic Riccati equation for the controller part, as given below:

 \begin{align} P =&\bar{A}^{\rm T}P\bar{A}+\bar{C}_{2}^{\rm T}\bar{C}_{2}-(\bar{C}_{2}^{\rm T}\bar{D} _{2}+\bar{A}^{\rm T}P\bar{B})\times \notag \\ &(\bar{D}_{2}^{\rm T}\bar{D}_{2}+\bar{B}^{\rm T}P\bar{B})^{-1}(\bar{D}_{2}^{\rm T}\bar{C }_{2}+\bar{B}^{\rm T}P\bar{A}). \label{AlgRic1} \end{align} (A9)

Defining

 \begin{align} {M}=\left[ \begin{array}{cc} \bar{A}^{\rm T}P\bar{A}-P+\bar{C}_2^{\rm T}\bar{C}_2&\bar{A}^{\rm T}P\bar{B}+\bar{C}_2^{\rm T}\bar{ D}_2 \\ \bar{B}^{\rm T}P\bar{A}+\bar{D}_2^{\rm T}\bar{C}_2&\bar{B}^{\rm T}P\bar{B}+\bar{D}_2^{\rm T}\bar{ D}_2 \end{array} \right] \end{align} (A10)

the stabilising controller $\bar{K}\in {\bf R}^{1\times n}$, from

 \begin{align} {M}\left[ \begin{array}{c} 0 \\ I \end{array} \right]=0 \end{align} (A11)

for a positive definite matrix $P\in {\bf R}^{n \times n}$, is given by

 \begin{align} \bar{K}=(\bar{B}^{\rm T}P\bar{B}+\bar{D}_2^{\rm T}\bar{D}_2)^{-1}(\bar{B}^{\rm T}P\bar{A}+ \bar{D}_2^{\rm T}\bar{C}_2). \end{align} (A12)

The property of duality is utilised for the estimation part of the ${H}_2$ optimal problem.

Table ${paragraph.serialNum} Accordingly, the algebraic Riccati equation for the estimation part is obtained from (A9) and the above given duality, as follows $ \begin{align} Q =&\bar{\bar{A}}Q\bar{\bar{A}}^{\rm T}+\bar{E}\bar{E}^{\rm T}-(\bar{\bar{A}}Q\bar{C% }_{1}^{\rm T}+\bar{E}\bar{D}_{1}^{\rm T})\times \notag \\ &(\bar{D}_{1}\bar{D}_{1}^{\rm T}+\bar{C}_{1}Q\bar{C}_{1}^{\rm T})^{-1}(\bar{D}_{1} \bar{E}^{\rm T}+\bar{C}_{1}Q\bar{\bar{A}}^{\rm T}) \end{align} $(A13) where$Q\in {\bf R}^{n\times n}$is a positive (semi-)definite matrix. And the corresponding estimator gain$\bar{L}\in {\bf R}^{n\times n}$is given as $ \begin{align} \bar{L}=(\bar{\bar{A}}Q\bar{C}_{1}^{\rm T}+\bar{E}\bar{D}_{1}^{\rm T})(\bar{D}_{1} \bar{D}_{1}^{\rm T}+\bar{C}_{1}Q\bar{C}_{1}^{\rm T})^{-1}. \end{align} $(A14) Appendix B Mechanistic model of a bioreactor The mechanistic model of a bioreactor is given by the following differential equations [14, 15] $ \begin{align} \frac{{\rm d}X}{{\rm d}t} =&-DX+\mu X \end{align} $(B1) $ \begin{align} \frac{{\rm d}S}{{\rm d}t} =&D(S_{f}-S)-\frac{1}{Y_{\frac{X}{S}}}\mu X \end{align} $(B2) $ \begin{align} \frac{{\rm d}P}{{\rm d}t} =&-DP+(\alpha \mu +\beta )X \end{align} $(B3) where$X$represents effluent cell-mass or biomass concentration,$S$represents substrate concentration and$P$denotes product concentration. Product concentration ($P$) and the cell-mass concentration ($X$) are measured process outputs, while dilution rate ($D$) and the feed substrate concentration ($S_{f}$) are the process inputs.$Y_{\frac{X}{S}}$represents the cell-mass yield and$\alpha $and$\beta $are the yield parameters for the product. The specific growth rate,$\mu $, is given by $ \begin{align} \mu =\frac{\mu _{m}(1-\dfrac{P}{P_{m}})}{K_{m}+S+\dfrac{S^{2}}{K_{i}}} \end{align} $(B4) where$\mu _{m}$represents maximum specific growth rate,$P_{m}$,$K_{m}$and$K_{i}$are the product saturation constant, substrate saturation constant and substrate inhibition constant, respectively. The nominal model parameters are given in the following table. Table${paragraph.serialNum}

References
 [1] Pakkiriswamy Aadaleesan and Prabirkumar Saha. A Nash game approach to mixed H2/H∞ model predictive control: Part 1-State feedback linear system. International Journal of Dynamics and Control, DOI: 10.1007/s40435-016-0261-y,to be published. [2] Pakkiriswamy Aadaleesan and Prabirkumar Saha. A Nash game approach to mixed H2/H∞ model predictive control: Part 2-Stability and Robustness. International Journal of Dynamics and Control, DOI: 10.1007/s40435-016-0259-5,to be pulbished. [3] Aadaleesan P., Miglani N., Sharma R., Saha P.. Nonlinear system identification using Wiener type LaguerreWavelet network model. Chemical Engineering Science, vol.63, no.15, pp.3932-3941, 2008. DOI:10.1016/j.ces.2008.04.043 [4] P. Aadaleesan. Nash game based mixed H2/H∞ model predictive control applied with Laguerre-Wavelet network model, Ph. D. dissertation, ⅡT Guwahati, India, 2011. [5] A. Pakkirisamy, P. Saha. Output feedback mixed H2/H∞ model predictive control: Observer based approach. In Proceedings of Annual Meeting of American Institute of Chemical Engineers, Salt Lake City, USA, 2010. https://www.researchgate.net/publication/267320581_Mixed_H2H_Model_Predictive_Control_for_Unstable_and_Non-Minimum_Constrained_Processes [6] Mayne D. Q., Raković S. V., Findeisen R., Allgöwer F.. Robust output feedback model predictive control of constrained linear systems. Automatica, vol.42, no.7, pp.1217-1222, 2006. DOI:10.1016/j.automatica.2006.03.005 [7] Trentelman H. L., Stoorvogel A. A.. Sampled-data and discrete-time H2 optimal control. SIAM Journal on Control and Optimization, vol.33, no.3, pp.834-862, 1995. DOI:10.1137/S0363012992241995 [8] Zhang W. H., Huang Y. L., Xie L. H.. Infinite horizon stochastic H2/H∞ control for discrete-time systems with state and disturbance dependent noise. Automatica, vol.44, no.9, pp.2306-2316, 2008. DOI:10.1016/j.automatica.2008.01.028 [9] Seki H., Ogawa M., Ooyama S., Akamatsu K., Ohshima M., Yang W.. Industrial application of a nonlinear model predictive control to polymerization reactors. Control Engineering Practice, vol.9, no.8, pp.819-828, 2001. DOI:10.1016/S0967-0661(01)00046-6 [10] Aadaleesan P., Saha P.. Nonlinear system identification using Laguerre wavelet models. Chemical Product and Process Modeling, vol.3, no.2, pp.Article number 3, 2008. [11] A. Griewank. Evaluating Derivatives: Principles and Techniques of Algorithmic Differentiation, Philadelphia, USA: Society for Industrial and Applied Mathematics, 2000. [12] S. A. Forth. User guide for MAD-A MATLAB Automatic Differentiation Toolbox, Applied Mathematics and Operational Research Report AMOR 2001/5, Cranfield University, Swindon, SN68LA, UK, June 2001. [13] Shampine L. F.. Accurate numerical derivatives in MATLAACM B. Transactions on Mathematical Software, vol.33, no.4, pp.Article number 26, 2007. DOI:10.1145/1268776 [14] Henson M. A., Seborg D. E.. Nonlinear control strategies for continuous fermenters. Chemical Engineering Science, vol.47, no.4, pp.821-835, 1992. DOI:10.1016/0009-2509(92)80270-M [15] Saha P., Krishnan S. H., Rao V. S. R., Patwardhan S. C.. Modeling and predictive control of MIMO nonlinear systems using Wiener-Laguerre models. Chemical Engineering Communications, vol.191, no.8, pp.1083-1119, 2004. DOI:10.1080/0098644049276452 [16] P. Saha. Nonlinear Model Predictive Control of Chemical Processes Using Laguerre Models, Ph. D. dissertation, ⅡT Madras, India, 1999. [17] Ding B., Huang B.. Output feedback model predictive control for nonlinear systems represented by hammersteinwiener model. IET Control Theory & Applications, vol.1, no.5, pp.1302-1310, 2007.