# Experimental Application of Predictive Controllers

Authors: C. H. F. Silva; H. M. Henrique; L. C. Oliveira-Lopes
Journal: Journal of Control Science and Engineering
 (2012)

---

## Abstract

Model predictive control (MPC) has been used successfully in industry. The basic characteristic of these algorithms is the formulation of an optimization problem in order to compute the sequence of control moves that minimize a performance function on the  time horizon with the best information available at each instant, taking into account operation and plant model constraints. The classical algorithms Infinite Horizon Model Predictive Control (IHMPC) and Model Predictive Control with Reference System (RSMPC) were used for the experimental application in the multivariable control of the pilot plant (level and pH). The simulations and experimental results indicate the applicability and limitation of the control technique.


---

## Body

## 1. Introduction

The process control is related to the application of the automatic control principles to industrial processes. The globalization effect upon industries brought a perception of the importance associated to the product quality over organization profits. Because of that, the process control has been more and more demanded and explored, in order not only to assure that the process is according to acceptable performance levels but also to address legal requirements in terms of safety and quality of the products [1]. In this context, predictive controllers are able to deal with system requirements in a proper way and simple to be implemented. One of the biggest concerns on control theory is related to the stability of closed loop of systems. It is natural that only stable closed-loop response is considered for real implementation. A representative algorithm of this controller class is Infinite Horizon Model Predictive Control (IHMPC) [2].During the closed-loop project, there are some theoretical tools that can be incorporated into the controllers and aggregate desirable characteristics. This is the case of predictive controllers that incorporate a reference path in their formulation (RSMPC) [3]. In this case, the value of control action is computed from a QDMC (Quadratic Dynamic Matrix Control) problem that has a first-order reference system incorporated directly into the formulation. Consequently, there is an effort for eliminating undesirable situations of too high speed and for balancing the system response.Most industrial chemical processes have a nonlinear characteristic, although linear controllers are used in many of these systems. The great advantage of this approximation is having an analytical solution to the control problem and also the low computational complexity. This kind of approach is very common in the industry practice.The pH system is used as a benchmark for applications in process control, mainly because of its strong nonlinear behavior. For the experimental application case addressed in this work, it is a control of multiple inputs and multiple outputs (MIMO), in which inputs are acid and base flows and outputs are reactor level and stream pH output. This system presents not only a nonlinear feature, but also an interaction between inputs and outputs linked to the system directionality. In this kind of system, classical controllers do not work to keep the stability and required performance. In addition, MPC solves the problem of open loop optimal control in real time, considering the current state of the system to be controlled as much as a policy feedback offline. As can be seen in the experimental data, the multivariable system could not be controlled with optimal classical controllers. That is why it was chosen to implement IHMPC and RSMPC predictive controllers in real time [4].Section2 presents the control algorithms used in the implementation described in this work. In Section 3 there is a description of the experimental system. Section 4 presents the simulations with the controllers and experimental results of the closed loop. Section 5 shows the main conclusions of this paper.


## 2. Model Predictive Control

The model predictive control (MPC) theory was originated in the end of 1970s decade and developed considerably ever since [5–10]. MPC became an important control strategy of industrial applications, mainly due to great success of its implementations in the petrochemical industry. The reason for MPC being so popular might be due to its ability to deal with quite difficult control challenges, although other characteristics such as its ability to control constrained multivariable plants, have been claimed as one of the most important features [6]. The ability to handle input, output, and internal state constraints can be assigned as the most significant contribution to the many successful industrial applications [11–13]. Other benefits to consider are the possibility of embedding safety limits, the possibility of obtaining advantages in using it in highly nonlinear plants as well as in time-varying plants, the consideration of process noise in the formulation, the high level of control performance, reduced maintenance requirements, and improved flexibility and agility [14].The various predictive control algorithms differ from each other depending on the way the predictive model is used to represent the process, the noise description, and on the cost function to be minimized. There are many control predictive applications well succeeded not only in chemical industry but also in other areas [15]. The use of state-space models, in spite of other formulations, was responsible for a substantial maturing of the predictive control theory during the 1990s decade [16]. The state-space formulation not only allows the application of linear system theorems already known, but also it makes easier to generalize for more complex cases. In this situation, the MPC controller can be understood as a compensator based on a state observer and its stability, performance, and robustness are determined by the observer poles, established directly by parameter adjustment, and by regulator poles, determined by performance horizons and weights [17]. Figure 1 shows the receding horizon strategy, which is the center of the MPC theory.Figure 1
MPC and the receding horizon strategy.Figure2, based on Richalet [18] and Richards and How [19], shows interfacial areas of the MPC control developments.Figure 2
MPC and related areas.Considering its application to MIMO processes, MPC deals directly with coordination and balancing interactions between inputs and outputs and also with associated constraints [20]. In spite of being really effective in suitable situations, MPC still has limitations as, for instance, operation difficulties, high maintenance cost (as it demands specialized work) and flexibility loss, which can result in a controller weakness. These limitations are not connected only to the algorithm itself, but also to the necessity of a plant model, for which is necessary maintenance. In practical terms, the limits of MPC applicability and performance should not be connected to the algorithm deficiencies, but to the issues linked to the modeling difficulties, sensor adequacy, and insufficient robustness in the presence of faults [15, 21].


### 2.1. IHMPC

The IHMPC formulation was introduced by Muske and Rawlings [2]. For simplicity the algorithm will be presented only for stable systems. For unstable systems see the reference section [2]. The discrete dynamical system used is shown in (1) in which y is the vector output, u is the vector input, and x is the states vector:(1)xk+1=Axk+1+Buk,yk=Cxk.
The receding horizon regulator is based on the minimization of the infinite horizon open-loop quadratic objective function at time k (2). Q is a symmetric positive semidefinite penalty matrix on the outputs. R is a symmetric positive semidefinite penalty matrix on the inputs. S is a symmetric positive semidefinite penalty matrix on the rate of the input change with Δu=uk+i-uk. The vector uNcontains the N future open-loop control moves (3). The infinite horizon open-loop objective function can be expressed as a finite horizon open-loop objective. For stable system, Q̅ is defined as the infinite sum (4). This infinite sum can be determined from the solution of the discrete Lyapunov equation (5). Using simple algebraic manipulation, the quadratic objective is show in (6). The matrices H, G, and F are showed in (7)–(9):(2)Jk=min⁡uN∑i=0∞yk+iTQyk+i+uk+iTRuk+i+Δuk+iTSΔuk+i,(3)uN=[ukuk+1⋮uk+N-1],(4)Q̅=∑i=0∞ATiCTQCAi,(5)Q̅=CTQC+ATQ̅A,(6)min⁡ΦkuN=min⁡uNuNTHuN+2uNT(Gxk-Fuk-1),(7)H=[BTQ̅B+R+2SBTATQ̅B-S⋯BTATN-1Q̅BBTQ̅AB-SBTQ̅B+R+2S⋯BTATN-2Q̅B⋮⋮⋱⋮BTQ̅AN-1BBTQ̅AN-2B⋯BTQ̅B+R+2S],(8)G=[BTQ̅ABTQ̅A2⋮BTQ̅AN],(9)F=[S0⋮0].


### 2.2. RSMPC

The RSMPC strategy is basically a transformation of an optimization problem with constraints into a quadratic programming problem [3]. The RSMPC development is presented next. Consider a general system described by (10). By linearizing (10) with the Taylor’s series expansion around the point immediately before the current operational point, one obtains (11). The matrices are defined in (12)–(16). Equation (11) is integrated from t to t+Δt and by assuming uk constant between samples, one can write (17):(10)dxdt=f(x,u),y=h(x),(11)dxdt=Ak-1x+Bk-1u+fk-1,y=Ck-1x+hk-1,(12)Ak-1=(∂f∂x)|x=xk-1,u=uk-1,(13)Bk-1=(∂f∂u)|x=xk-1,u=uk-1,(14)Ck-1=(∂h∂x)|x=xk-1,(15)fk-1=f(x,u)|x=xk-1,u=uk-1,(16)hk-1=h|x=xk-1,(17)x̂k+1=Φxk+ΨBuk+Ψfk-1,ŷk=Cxk+hk-1,
where
(18)Φ=eAΔt,Ψ=A-1(eAΔt-I)B,Ω=A-1(eAΔt-I).
Since A is not singular, then (17) can be written for each prediction instant from k=1 to k=P, where P is the prediction horizon and N is the control horizon, with P≥N and Δuk+j=0 to N≤j≤P. The resulting set of equation is showed below, where(19)ŷ=ΓΔu+γ,(20)ŷ=[ŷk+1ŷk+2⋯ŷk+N⋯ŷk+P],(21)Δu=[Δuk+1Δuk+2⋯Δuk+N].
The controller is designed in order to transform the closed-loop behavior in a fist-order system behavior according to (24). To overcome the problem of infeasible solutions of the MPC optimization problem, a slack variable (λ) is introduced in order to allow the system to deviate from reference system and to satisfy hard constraints [3] resulting in (25). This equation can be written for each prediction instant (26): (22)Γ=[CΨ0⋯0⋮⋮⋮C(Φ+I)ΨCΨ⋯0⋮⋮⋮C(∑i=1NΦi-1)ΨC(∑i=1N-1Φi-1)Ψ⋯CΨ⋮⋮⋮C(∑i=1PΦi-1)ΨC(∑i=1P-1Φi-1)Ψ⋯C(∑i=1P-NΦi-1)Ψ],(23)γ=[yk+CΦΔxkyk+C(∑i=12Φi)Δxk⋮yk+C(∑i=1NΦi)Δxk⋮yk+C(∑i=1PΦi)Δxk],(24)dydt=K(ySP-y),(25)Cfk-1+CAk-1xk+CBk-1uk+λk=K(ySP-y),(26)Dυ=b,
where
(27)υ=[Δuk⋯Δuk+Nλk⋯λk+P],(28)b=[KykSP-(CA+KC)xk-Cfk-1Kyk+1SP-(CA+KC)ΦxkKyk+2SP-(CA+KC)Φ2xk⋮Kyk+PSP-(CA+KC)ΦPxk],(29)D=[CB0⋯0I000⋯0(CA+KC)ΨCB⋯0-II00⋯0(CA+KC)ΨC(CA+KC)Ψ⋯00-II0⋯0⋮⋮⋱⋮⋮⋮⋮⋮⋱⋮(CA+KC)ΦP-1Ψ(CA+KC)ΦP-2Ψ⋯(CA+KC)ΦP-M-1Ψ0000⋯I].The cost function can be written as (30). This equation can also be reorganized as a quadratic programming problem (31):(30)min⁡Δu(k),…,Δu(k+Nm)λ(k),…,λ(k+Np)J=12(ΔuTRΔu+λTSλ),(31)min⁡υ(k),…,υ(k+Nm+Np+2)⁡J=12υTευsubjectto:Dυ=bumin⁡k+j≤uk+j≤umax⁡k+j-|Δumin⁡|k+j≤Δuk+j≤|Δumax⁡|k+jj=0,…,P
in which(32)ε=[R[0][0]S].
The controller cannot eliminate offset unless (33) is handled [3]:(33)Cfk-1≈yk-yk-1Δt.


## 2.1. IHMPC

The IHMPC formulation was introduced by Muske and Rawlings [2]. For simplicity the algorithm will be presented only for stable systems. For unstable systems see the reference section [2]. The discrete dynamical system used is shown in (1) in which y is the vector output, u is the vector input, and x is the states vector:(1)xk+1=Axk+1+Buk,yk=Cxk.
The receding horizon regulator is based on the minimization of the infinite horizon open-loop quadratic objective function at time k (2). Q is a symmetric positive semidefinite penalty matrix on the outputs. R is a symmetric positive semidefinite penalty matrix on the inputs. S is a symmetric positive semidefinite penalty matrix on the rate of the input change with Δu=uk+i-uk. The vector uNcontains the N future open-loop control moves (3). The infinite horizon open-loop objective function can be expressed as a finite horizon open-loop objective. For stable system, Q̅ is defined as the infinite sum (4). This infinite sum can be determined from the solution of the discrete Lyapunov equation (5). Using simple algebraic manipulation, the quadratic objective is show in (6). The matrices H, G, and F are showed in (7)–(9):(2)Jk=min⁡uN∑i=0∞yk+iTQyk+i+uk+iTRuk+i+Δuk+iTSΔuk+i,(3)uN=[ukuk+1⋮uk+N-1],(4)Q̅=∑i=0∞ATiCTQCAi,(5)Q̅=CTQC+ATQ̅A,(6)min⁡ΦkuN=min⁡uNuNTHuN+2uNT(Gxk-Fuk-1),(7)H=[BTQ̅B+R+2SBTATQ̅B-S⋯BTATN-1Q̅BBTQ̅AB-SBTQ̅B+R+2S⋯BTATN-2Q̅B⋮⋮⋱⋮BTQ̅AN-1BBTQ̅AN-2B⋯BTQ̅B+R+2S],(8)G=[BTQ̅ABTQ̅A2⋮BTQ̅AN],(9)F=[S0⋮0].


## 2.2. RSMPC

The RSMPC strategy is basically a transformation of an optimization problem with constraints into a quadratic programming problem [3]. The RSMPC development is presented next. Consider a general system described by (10). By linearizing (10) with the Taylor’s series expansion around the point immediately before the current operational point, one obtains (11). The matrices are defined in (12)–(16). Equation (11) is integrated from t to t+Δt and by assuming uk constant between samples, one can write (17):(10)dxdt=f(x,u),y=h(x),(11)dxdt=Ak-1x+Bk-1u+fk-1,y=Ck-1x+hk-1,(12)Ak-1=(∂f∂x)|x=xk-1,u=uk-1,(13)Bk-1=(∂f∂u)|x=xk-1,u=uk-1,(14)Ck-1=(∂h∂x)|x=xk-1,(15)fk-1=f(x,u)|x=xk-1,u=uk-1,(16)hk-1=h|x=xk-1,(17)x̂k+1=Φxk+ΨBuk+Ψfk-1,ŷk=Cxk+hk-1,
where
(18)Φ=eAΔt,Ψ=A-1(eAΔt-I)B,Ω=A-1(eAΔt-I).
Since A is not singular, then (17) can be written for each prediction instant from k=1 to k=P, where P is the prediction horizon and N is the control horizon, with P≥N and Δuk+j=0 to N≤j≤P. The resulting set of equation is showed below, where(19)ŷ=ΓΔu+γ,(20)ŷ=[ŷk+1ŷk+2⋯ŷk+N⋯ŷk+P],(21)Δu=[Δuk+1Δuk+2⋯Δuk+N].
The controller is designed in order to transform the closed-loop behavior in a fist-order system behavior according to (24). To overcome the problem of infeasible solutions of the MPC optimization problem, a slack variable (λ) is introduced in order to allow the system to deviate from reference system and to satisfy hard constraints [3] resulting in (25). This equation can be written for each prediction instant (26): (22)Γ=[CΨ0⋯0⋮⋮⋮C(Φ+I)ΨCΨ⋯0⋮⋮⋮C(∑i=1NΦi-1)ΨC(∑i=1N-1Φi-1)Ψ⋯CΨ⋮⋮⋮C(∑i=1PΦi-1)ΨC(∑i=1P-1Φi-1)Ψ⋯C(∑i=1P-NΦi-1)Ψ],(23)γ=[yk+CΦΔxkyk+C(∑i=12Φi)Δxk⋮yk+C(∑i=1NΦi)Δxk⋮yk+C(∑i=1PΦi)Δxk],(24)dydt=K(ySP-y),(25)Cfk-1+CAk-1xk+CBk-1uk+λk=K(ySP-y),(26)Dυ=b,
where
(27)υ=[Δuk⋯Δuk+Nλk⋯λk+P],(28)b=[KykSP-(CA+KC)xk-Cfk-1Kyk+1SP-(CA+KC)ΦxkKyk+2SP-(CA+KC)Φ2xk⋮Kyk+PSP-(CA+KC)ΦPxk],(29)D=[CB0⋯0I000⋯0(CA+KC)ΨCB⋯0-II00⋯0(CA+KC)ΨC(CA+KC)Ψ⋯00-II0⋯0⋮⋮⋱⋮⋮⋮⋮⋮⋱⋮(CA+KC)ΦP-1Ψ(CA+KC)ΦP-2Ψ⋯(CA+KC)ΦP-M-1Ψ0000⋯I].The cost function can be written as (30). This equation can also be reorganized as a quadratic programming problem (31):(30)min⁡Δu(k),…,Δu(k+Nm)λ(k),…,λ(k+Np)J=12(ΔuTRΔu+λTSλ),(31)min⁡υ(k),…,υ(k+Nm+Np+2)⁡J=12υTευsubjectto:Dυ=bumin⁡k+j≤uk+j≤umax⁡k+j-|Δumin⁡|k+j≤Δuk+j≤|Δumax⁡|k+jj=0,…,P
in which(32)ε=[R[0][0]S].
The controller cannot eliminate offset unless (33) is handled [3]:(33)Cfk-1≈yk-yk-1Δt.


## 3. The Experimental System

The experimental system built was a neutralization process that occurs in a shaken reactor tank. Figure3 presents a photo of the experimental system.Figure 3
Experimental system.Figure4 shows a diagram indicating how the system works. The numbers indicate acid flow (1), base flow (2), buffer flow (3), and output flow (4). This number is used in the indexes of the modeling of the experimental system. The flows of acid (q1) and base (q2) are the manipulated variables (u); the flow of buffer solution (q3) is the unmeasured or measured disturbance (d), depending on how it is treated in the control formulation. For the case of this study, the buffer solution flow was used only to make simulation tests. The level height of the reactor (h) and the output flow pH (q4) are the controlled variables or the outputs (y).Figure 4
Diagram representation of the experimental system.The system was specially chosen for this study because of its strong nonlinear dynamic behavior and low operational cost, considering the consumption of inputs as electrical energy and raw material [3, 22]. The system becomes interesting under the control perspective, as it is a multivariable and nonlinear process in which the static gain suffers strong variations inside the flow operation band. Figure 5 presents the system titration curve. It can be observed that the process pH changes from 4 to 9 in the [0.87;1.22] interval of the ratio q2/q1. This change in the process gain is just what makes the neutralization problem difficult to solve under the control perspective, as small modifications on the relation q2/q1 can cause big or small pH variations, depending on where the operating point is in the operational region. It has to be emphasized that, in an analytical process of titration, only one drop is enough to change the pH from 4 to 9. Thus, it can be inferred that the sensitivity of the dynamic process associated to this offers a challenge problem for control in real time. It can be noticed that, once a buffer solution is added, the nonlinear behavior of the system smoothes up, so that the controller action gets easier, according to what the simulation curves with this solution indicate. If the system was linear, the relation given by pH ×q2/q1 plots would be a line whose approximation can be realized by the simulated curve of larger amount of buffer solution. Montandon [3] presented results for the closed-loop system using pH control of the acid flow with classical PID controller. The conclusion was that this implementation was not able to reject disturbances in a satisfactory manner. During the preparation of this work it has been studied certain PID testing for the MIMO case, but the answers proved to be mostly unstable and hard to tune.Figure 5
Titration Curves of the neutralization system.


### 3.1. Phenomenological Model

Hall [23] developed the physical model of this process, which is based on the hypotheses of the perfect mixture, constant density, and total solubility of the present ions. The chemical reactions involved in the acid-base neutralization (HNO3–NaOH) are shown in (34)–(38). In order to have a complete model it will be considered the presence of a buffer (NaHCO3). However, for the experimental results, this solution was not used:(34)HNO3⟶H++NO3-,(35)NaOH⟶Na++OH-,(36)H2O⟷kwH++OH-,(37)NaHCO3+H2O⟶NaOH+H2CO3,H2CO3⟷ka1H++HCO3-,(38)HCO3-⟷ka2H++CO3-2,
where(39)ka1=⌊HCO3-⌋⌊H+⌋[H2CO3],(40)ka2=⌊CO3-⌋⌊H+⌋[HCO3-].
The amounts of Wa and Wb are called as invariants, because their concentrations are not affected along the reaction. The reactions are fast enough to allow the system to be considered as in equilibrium. Then, the equilibrium reactions can be used to determine the concentration of the hydrogen (H) ions through the concentration of the reaction invariants. Equation (41) gives the equilibrium concentration. According to Gustafsson and Waller [24], there are two reaction invariants defined for the ith stream and presented in (42) and (43), which, when are combined, result in an implicit relation of H, Wa, and Wb, as (44) shows:(41)kw=⌊H+⌋⌊OH-⌋,(42)Wai=⌊H+⌋i-⌊OH-⌋i-⌊HCO3-⌋i-2⌊CO3-2⌋i,(43)Wbi=[H2CO3]i+[HCO3-]i-2[CO3-2]i,(44)Wa=H-kwH-Wb(ka1/H)+2(ka1ka2)/H21+(ka1/H)+(ka1ka2/H2).
Making the mass balance in the reactor, together with the invariant equations and considering that the density is constant by hypothesis, will result in a differential equation system (45)–(47):(45)dhdt=1Ar(q1+q2+q3-q4),(46)dWadt=1Vr(q1(Wa1-Wa)+q2(Wa2-Wa)+q3(Wa3-Wa)),(47)dWbdt=1Vr(q1(Wb1-Wb)+q2(Wb2-Wb),+q3(Wb3-Wb)),
where Ar is the reactor area, Vr is the reaction volume; h is the solution height in the reactor; Wa is the acid reaction invariant; Wb is the base reaction invariant.The flow output of the system is driven by gravity. The output flow(q4) is done by a globe-type valve (48). This equation shows the relation between output flow and the reactor height and the parameters to be estimated (cv and p7). The valve parameters were estimated by using the model and experimental data from open loop. After this process of defining the parameters of the valve it was locked so as to avoid modifying these parameters. Nominal operation conditions are presented in Table 1:(48)q4=cvhp7.Table 1
System nominal parameters.
VariablesSymbolNominal valuesLevelH25 cmAreaAr168.38 cm2VolumeVr4209.67 cm3Acid flowq111.9130 mL/sBase flowq211.8235 mL/sBuffer solution flowq30.01 mL/spHpH7.0Acid conc. inq1[HNO3]3.510×10-03 MBase conc. inq2[NaOH]3.510×10-03 MFor estimation and validation of the model, the plant operation must be done through the inputs that excite all the dynamic modes and in a large frequency band, so that the outputs have enough information for that procedure. These are the most usual excitation types: ramp, step, and impulse. However, in order to gather information in a broader range of frequency, random inputs made of sequences of steps with variable amplitude and duration are applied, as it is indicated by Montandon [3]. The sampling time for the process was 10 s. The acid and base flows vary as q1∈[9.85613.310], q2∈[9.97713.215], with a step probability equal to 0,8, in which the duration of each step was 120 s. The flows in which the variations occurred were 12 mL/s each, in order to keep the flow sum inside the limit given by q1+q2∈[21.524.5]. The values presented were defined in such a way that, during the operation, the reactor did not get empty nor presents overflow. A set of experimental data composed by 6000 acquisition points was used. The estimation procedure used half of the data and the validation employed the remaining data, in order to assure that the model would have the ability for predicting unknown data. The estimation is a nonlinear solution of the approximation curve that uses the least square method with a first trial given by a random value and a step equal to 15. The estimated values that resulted were cv=20.4477 and p7=0.0523. Figure 6 shows the error in predicting h. Figure 7 presents the set of input random data. Figure 8 presents the results of the simulation and the experimental ones for the closed-loop run, in which one can notice that the model is very representative of the process, with difficulties to predict around pH 7; this phenomenon was already expected, taking into account the titration curve of the system. Figure 9 shows the validation of the model by the prediction of the output that comes very close to the experimental data.Figure 6
Error in predicting the level of the reactor.Figure 7
Random inputs for parameter estimation.Figure 8
Open-loop run behavior (Mod indicates that it is a result of model simulation).Figure 9
Real datax simulated data (Mod).


### 3.2. State-Space Model

The state-space approach resulting from the algebraic manipulation is given by (49)–(55). The terms dWadz and dWbdz are defined in (56) and (57). Equation (58) presents dydz relation:(49)x=[hWaWb],(50)u=[q1q2],(51)y=[hpH],(52)d=[q3],(53)A=[-cvp7hp7-100-q1(Wa1-Wa)+q2(Wa2-Wa)+q3(Wa3-Wa)Arh2-q1-q2-q3Arh0-q1(Wb1-Wb)+q2(Wb2-Wb)+q3(Wb3-Wb)Arh20-q1-q2-q3Arh],(54)B=[1Ar1Ar(Wa1-Wa)Arh(Wa2-Wa)Arh(Wb1-Wb)Arh(Wb2-Wb)Arh],G=[1Ar(Wa3-Wa)Arh(Wb3-Wb)Arh],(55)C=[1000dydzdWadzdydzdWdz],(56)dWadz=1+kwH2-Wb(-(ka1/H2)-(4ka1⋅ka2/H3))1+(ka1/H)+(ka1⋅ka2/H2)+Wb((ka1/H)+(2⋅ka1.ka2/H2))(-(ka1/H2)+(2⋅ka1⋅ka2/H3))(1+(ka1/H)+(ka1⋅ka2/H2))2,(57)dWdz=(1+(kw/H2))(1+(ka1/H)+(ka1⋅ka2/H2))(ka1/H)+(2⋅ka1⋅ka2/H2)-(H-(kw/H)-Wa)(1+(ka1/H)+(ka1⋅ka2/H2))(-(ka1/H2)-(4⋅ka1⋅ka2/H3))((ka1/H)+(2⋅ka1⋅ka2/H2))2+(H-(kw/H)-Wa)(-(ka1/H2)-(2⋅ka1⋅ka2/H3))(ka1/H)+(2⋅ka1⋅ka2/H2),(58)dydz=-log10eH.


## 3.1. Phenomenological Model

Hall [23] developed the physical model of this process, which is based on the hypotheses of the perfect mixture, constant density, and total solubility of the present ions. The chemical reactions involved in the acid-base neutralization (HNO3–NaOH) are shown in (34)–(38). In order to have a complete model it will be considered the presence of a buffer (NaHCO3). However, for the experimental results, this solution was not used:(34)HNO3⟶H++NO3-,(35)NaOH⟶Na++OH-,(36)H2O⟷kwH++OH-,(37)NaHCO3+H2O⟶NaOH+H2CO3,H2CO3⟷ka1H++HCO3-,(38)HCO3-⟷ka2H++CO3-2,
where(39)ka1=⌊HCO3-⌋⌊H+⌋[H2CO3],(40)ka2=⌊CO3-⌋⌊H+⌋[HCO3-].
The amounts of Wa and Wb are called as invariants, because their concentrations are not affected along the reaction. The reactions are fast enough to allow the system to be considered as in equilibrium. Then, the equilibrium reactions can be used to determine the concentration of the hydrogen (H) ions through the concentration of the reaction invariants. Equation (41) gives the equilibrium concentration. According to Gustafsson and Waller [24], there are two reaction invariants defined for the ith stream and presented in (42) and (43), which, when are combined, result in an implicit relation of H, Wa, and Wb, as (44) shows:(41)kw=⌊H+⌋⌊OH-⌋,(42)Wai=⌊H+⌋i-⌊OH-⌋i-⌊HCO3-⌋i-2⌊CO3-2⌋i,(43)Wbi=[H2CO3]i+[HCO3-]i-2[CO3-2]i,(44)Wa=H-kwH-Wb(ka1/H)+2(ka1ka2)/H21+(ka1/H)+(ka1ka2/H2).
Making the mass balance in the reactor, together with the invariant equations and considering that the density is constant by hypothesis, will result in a differential equation system (45)–(47):(45)dhdt=1Ar(q1+q2+q3-q4),(46)dWadt=1Vr(q1(Wa1-Wa)+q2(Wa2-Wa)+q3(Wa3-Wa)),(47)dWbdt=1Vr(q1(Wb1-Wb)+q2(Wb2-Wb),+q3(Wb3-Wb)),
where Ar is the reactor area, Vr is the reaction volume; h is the solution height in the reactor; Wa is the acid reaction invariant; Wb is the base reaction invariant.The flow output of the system is driven by gravity. The output flow(q4) is done by a globe-type valve (48). This equation shows the relation between output flow and the reactor height and the parameters to be estimated (cv and p7). The valve parameters were estimated by using the model and experimental data from open loop. After this process of defining the parameters of the valve it was locked so as to avoid modifying these parameters. Nominal operation conditions are presented in Table 1:(48)q4=cvhp7.Table 1
System nominal parameters.
VariablesSymbolNominal valuesLevelH25 cmAreaAr168.38 cm2VolumeVr4209.67 cm3Acid flowq111.9130 mL/sBase flowq211.8235 mL/sBuffer solution flowq30.01 mL/spHpH7.0Acid conc. inq1[HNO3]3.510×10-03 MBase conc. inq2[NaOH]3.510×10-03 MFor estimation and validation of the model, the plant operation must be done through the inputs that excite all the dynamic modes and in a large frequency band, so that the outputs have enough information for that procedure. These are the most usual excitation types: ramp, step, and impulse. However, in order to gather information in a broader range of frequency, random inputs made of sequences of steps with variable amplitude and duration are applied, as it is indicated by Montandon [3]. The sampling time for the process was 10 s. The acid and base flows vary as q1∈[9.85613.310], q2∈[9.97713.215], with a step probability equal to 0,8, in which the duration of each step was 120 s. The flows in which the variations occurred were 12 mL/s each, in order to keep the flow sum inside the limit given by q1+q2∈[21.524.5]. The values presented were defined in such a way that, during the operation, the reactor did not get empty nor presents overflow. A set of experimental data composed by 6000 acquisition points was used. The estimation procedure used half of the data and the validation employed the remaining data, in order to assure that the model would have the ability for predicting unknown data. The estimation is a nonlinear solution of the approximation curve that uses the least square method with a first trial given by a random value and a step equal to 15. The estimated values that resulted were cv=20.4477 and p7=0.0523. Figure 6 shows the error in predicting h. Figure 7 presents the set of input random data. Figure 8 presents the results of the simulation and the experimental ones for the closed-loop run, in which one can notice that the model is very representative of the process, with difficulties to predict around pH 7; this phenomenon was already expected, taking into account the titration curve of the system. Figure 9 shows the validation of the model by the prediction of the output that comes very close to the experimental data.Figure 6
Error in predicting the level of the reactor.Figure 7
Random inputs for parameter estimation.Figure 8
Open-loop run behavior (Mod indicates that it is a result of model simulation).Figure 9
Real datax simulated data (Mod).


## 3.2. State-Space Model

The state-space approach resulting from the algebraic manipulation is given by (49)–(55). The terms dWadz and dWbdz are defined in (56) and (57). Equation (58) presents dydz relation:(49)x=[hWaWb],(50)u=[q1q2],(51)y=[hpH],(52)d=[q3],(53)A=[-cvp7hp7-100-q1(Wa1-Wa)+q2(Wa2-Wa)+q3(Wa3-Wa)Arh2-q1-q2-q3Arh0-q1(Wb1-Wb)+q2(Wb2-Wb)+q3(Wb3-Wb)Arh20-q1-q2-q3Arh],(54)B=[1Ar1Ar(Wa1-Wa)Arh(Wa2-Wa)Arh(Wb1-Wb)Arh(Wb2-Wb)Arh],G=[1Ar(Wa3-Wa)Arh(Wb3-Wb)Arh],(55)C=[1000dydzdWadzdydzdWdz],(56)dWadz=1+kwH2-Wb(-(ka1/H2)-(4ka1⋅ka2/H3))1+(ka1/H)+(ka1⋅ka2/H2)+Wb((ka1/H)+(2⋅ka1.ka2/H2))(-(ka1/H2)+(2⋅ka1⋅ka2/H3))(1+(ka1/H)+(ka1⋅ka2/H2))2,(57)dWdz=(1+(kw/H2))(1+(ka1/H)+(ka1⋅ka2/H2))(ka1/H)+(2⋅ka1⋅ka2/H2)-(H-(kw/H)-Wa)(1+(ka1/H)+(ka1⋅ka2/H2))(-(ka1/H2)-(4⋅ka1⋅ka2/H3))((ka1/H)+(2⋅ka1⋅ka2/H2))2+(H-(kw/H)-Wa)(-(ka1/H2)-(2⋅ka1⋅ka2/H3))(ka1/H)+(2⋅ka1⋅ka2/H2),(58)dydz=-log10eH.


## 4. Results and Discussions

The real-time implementation was done using the successive linearization around the operating point. In this case, the linear form results (11).A controller designed around a stationary state would not be able to control the system, because, for some regions, the stationary state would be very far from the desired point and so the plant model mismatch. Even with the successive linearization, the system control would not be successful, mainly in regions with high gain. There are several ways to deal with this problem. One of them is to expand the states, creating a new state matrix, and then the original structure of the optimization problem can be kept (59). The term related to the output linearization does originate problems for the controller, because it only maps state output and does not interfere in the systems dynamics as the term fk-1 does. Once this modification is done, it is possible to work on the controller previously defined. This option was chosen for real-time operation. Another form can be found in Reis [25]:(59)dxdt=[Ak-100fk-1][x1]+Bk-1u.
For the real application in the experimental system, which has nonlinear characteristics, it is necessary to perform successive linearization of the process model in order to apply linear algorithms. The dynamic evolution of the process results in a time-varying system (LTV). The equation of the IHMPC and RSMPC in an LTV form can be easily extrapolated from the formulation presented herein, and for simplicity its description is omitted.In the experimental system, there are the statesh, Wa, and Wb. For this case, the last two states are not measured, as it was shown in the modeling section of the system, so you need to estimate them (x̂). To this end, we used a Kalman filter [26], (60). For the IHMPC controller, we also used an open-loop observer (61), informing the model that the error (e) enters into the plant and model in order to minimize the interference of state estimation in the system behavior [27]:(60)x̂k+1∣k=Ax̂k∣k-1+Buk+L(yk-Cx̂k∣k-1),(61)e=yPlant-yModel.
RSMPC and IHMPC were implemented in real-time operation of the experimental plant, using the interface through LabVIEW and the controller computation, using the routines implemented in Matlab. In the literature there are several techniques of controller tuning. In this work Henson and Seborg [11] and Montadon’s [3] indications were adopted. Besides, a field refinement was done during the initial run in closed loop with each controller, making small adjustments in parameters, so that the system would deliver a better experimental performance in closed loop. Such results were omitted in this paper. Table 2 shows the parameters of each controller resulting from the simulation and used in the experimental run. The controllers were at first tuned with identity matrices and a control horizon equal to 10; besides, a tuning was done until achieving a satisfactory response. Figure 10 presents the structure of the experimental plant, and Figures 11 and 12 show the LabVIEW implementation.Table 2
Controller tuning.
ControllerParameterIHMPCN=10R=10I→100IQ=IS=10IΔumax=0.25 mL/sRSMPCK=0.005IR=I→R=100IS=10000IN=2P=10Δumax=0.5 mL/sGeneral parametersΔt=10 sumax=38 mL/sumin=4 mL/sFigure 10
Structure of operation—experimental plant.Figure 11
Frontal panel—IHMPC.Figure 12
Block diagram—IHMPC.Equation (62) shows that the calibration curves of the instruments. V0,V1,V2,V3, and V4 are channel voltages. The regression coefficients are equal to 1 [4, 28]:(62)q1=5.2788V0-9.4776,q2=5.0997V1-9.5541,q3=0.2090V2-0.3845,pH=1.8982V3-3.5583,h=5.6347V4-9.0838.
Setpoint variations were done to evaluate the ability of the controllers to lead with transitions. Starting the system with a pH equal to 7 and a reactor height equal to 18 cm, variations were done according to the vector in (63), considering that each step took 2000 seconds. For all controller simulations a random noise of average equal to zero was added with about 10% of the outputs (h and pH):(63)[hpHd]=[1818181818211815187107477777000000000].Figures13, 14, 15, and 16 present the closed-loop simulation that uses RSMPC. The performance shown is quite reasonable.Figure 13
Closed-loop simulation: RSMPC—pH (setpoint: SP and reference system: Ref).Figure 14
Closed-loop simulation: RSMPC—h (setpoint: SP and reference system: Ref).Figure 15
Closed-loop simulation: RSMPC—control moves.Figure 16
Closed-loop simulation: RSMPC—slack variables.The experimental run was done in two steps to preserve and manipulate the data separately. In the first part was made a variation in the pH setpoint and kept the level constant. In the second part it was done otherwise. Figures17, 18, 19, and 20 show the experimental run of the pH setpoint changes, and Figures 21, 22, 23, and 24 show the experimental run of level setpoint changes. The system was controlled in all runs. During the experiment, there was a need to increase the weight matrix of the control action (as indicated in Table 2); as the output response started to oscillate, it was not able to reach the desired setpoint. However, the overall response can be considered satisfactory.Figure 17
Experimental response RSMPC—pH (setpoint: SP and reference system: Ref).Figure 18
Experimental response: RSMPC—h (setpoint: SP and reference system: Ref).Figure 19
Experimental response: inputs.Figure 20
Experimental response: RSMPC—slack variables.Figure 21
Experimental response RSMPC—pH (setpoint: SP and reference system: Ref).Figure 22
Experimental response: RSMPC—h (setpoint: SP and reference system: Ref).Figure 23
Experimental response: inputs.Figure 24
Experimental response: RSMPC—slack variables.The results for the IHMPC simulation are presented in Figures25, 26, 27, and 28. The results are suitable to apply on line the control to the real process.Figure 25
Closed-loop simulation: IHMPC-pH.Figure 26
Closed-loop simulation: IHMPC-h.Figure 27
Closed-loop simulation: IHMPC-inputs.Figure 28
Closed-loop simulation: IHMPC-optimization time.Like RSMPC, experimental runs with IHMPC were done in two steps. Figures29, 30, and 31 show the experimental run of the pH setpoint changes, and Figures 32, 33, and 34 show the experimental run of level setpoint changes. This controller showed a performance superior to the RSMPC controller performance and with satisfactory and adequate response and performance.Figure 29
Experimental response: IHMPC-pH.Figure 30
Experimental response: IHMPC-h.Figure 31
Experimental response: IHMPC-inputs.Figure 32
Experimental response: IHMPC-pH.Figure 33
Experimental response: IHMPC-h.Figure 34
Experimental response: IHMPC-inputs.


## 5. Conclusions

Through this research it was able to better understand the use of predictive controllers in a real-time application, paving the way for research in the area. The experimental application led to an approximation of reality and industrial practice, experiencing some of the common problems and issues in the implementation of controllers. It was possible to verify the need for a mastery of techniques and concepts of process control, system modeling, parameter identification, scheduling, optimization, and, in addition, common sense engineering to solve the experimental problems.The control of such class of nonlinear processes is a very challenging area with many possibilities for development and that undoubtedly has importance and influence on the performance of process and in consequence results in their organizations.This work carried out the experimental application of two predictive controllers: IHMPC and RSMPC, both in simulation environment and in the experimental plant dealing with the level and pH control.The experimental plant was modeled and has the required parameters identified through the run in open loop. The simulation results were satisfactory and indicated that the model is representative of the real process and suitable for the control purposes that were aimed. The simulation answers of the system outputs subjected to the controllers were adequate and satisfactory. The theoretical and experimental responses of the IHMPC runs were satisfactory.The controller adjustment for real-time operations was suitable and feasible. The interference of matters related to leaking, process noises, and noises from electrical source and other problems associated to real-time application brought additional difficulties that demanded process knowledge of the process.The experimental application of robust controllers associated or not to control systems tolerant to failures and the use of online identification employing neural networks will be presented elsewhere.



---

*Source: 159072-2012-02-19.xml*