# Model Predictive Control with a Relaxed Cost Function for Constrained Linear Systems

Authors: David Sotelo; Antonio Favela-Contreras; Viacheslav V. Kalashnikov; Carlos Sotelo
Journal: Mathematical Problems in Engineering
 (2020)

---

## Abstract

The Model Predictive Control technique is widely used for optimizing the performance of constrained multi-input multi-output processes. However, due to its mathematical complexity and heavy computation effort, it is mainly suitable in processes with slow dynamics. Based on the Exact Penalization Theorem, this paper presents a discrete-time state-space Model Predictive Control strategy with a relaxed performance index, where the constraints are implicitly defined in the weighting matrices, computed at each sampling time. The performance validation for the Model Predictive Control strategy with the proposed relaxed cost function uses the simulation of a tape transport system and a jet transport aircraft during cruise flight. Without affecting the tracking performance, numerical results show that the execution time is notably decreased compared with two well-known discrete-time state-space Model Predictive Control strategies. This makes the proposed Model Predictive Control mainly suitable for constrained multivariable processes with fast dynamics.


---

## Body

## 1. Introduction

Model Predictive Control (MPC) for linear systems is now a well-established discipline providing stability, feasibility, and robustness [1–6]. Due to its inherent ability to take into account constraints and deal with multi-input multi-output variables [7–10], it has been applied in a wide range of applications, including chemical processes, industrial systems, energy, health, environment, and aerospace [11–16]. In [17], a robust MPC strategy is presented to handle the trajectory tracking problem for an underactuated two-wheeled inverted pendulum vehicle. Moreover, based on an MPC scheme, in [18], a control strategy is designed to an unmanned aerial vehicle for its automatic carrier landing system. Nevertheless, the computation complexity makes the multivariable MPC ineffectual for high speed applications where the controller must execute in a few milliseconds [19–23]. Moreover, the problem becomes much more complicated solving such an online constrained optimization problem by computing a numerical solver [24–26]. Several MPC techniques are used to overcome these problems. For instance, in [27], an explicit model predictive control moves major part of computation offline, which makes it enable to be implemented in real time for wide range of fast systems. Also, in [28], to reduce the online computational time, all the state trajectories are included in the optimal control problem as the constraints in the prediction horizon, then only a quadratic programming problem is solved. In [29], based on a mixed integer quadratic programming problem, the control input is calculated at each discrete time. In contrast to common MPC approaches, where an optimization toolbox is required, this work presents a relaxed performance index, in which the weighting matrices are computed online using the concept of Taylor series expansion and standard inverse distance weighting (IDW) functions. Then, tracking performance under input-output constraints is well obtained, lighter computation load is achieved, and execution time to solve a Quadratic Program (QP) is reduced. Thus, a computationally efficient constrained MPC for discrete-time state-space multivariable systems is obtained.The paper is organized as follows. Section2 gives the preliminaries of the proposed MPC strategy. Section 3 describes the proposed relaxed cost function. Section 4 presents a tape transport system and a jet transport aircraft as study cases. Simulation results show the performance of the proposed MPC strategy and the execution time improvement compared with two well-known MPC strategies. Finally, Section 5 discusses the conclusions. Acknowledgments and the list of references finish the paper.


## 2. Model Predictive Control Based on Discrete-Time State Space Model

This section presents a brief review of MPC based on discrete-time state-space model. The original controller is proposed by Alamir in [7]. In this previous work, considering the predictions of the states, the control action is obtained through the solution of a constrained optimization problem by using a cost function with constant weighting matrices. At each sampling time, an optimal control problem is solved whose results are computationally expensive. The system dynamics is denoted by the Linear Time Invariant (LTI) State-Space Model taking the following structure:(1)xk+1=Axk+Buk,yrk=Crxk,where x∈ℛn is the state vector, u∈ℛnu is the controlled input vector, yr∈ℛnr is the output vector, A∈ℛn×n is the state matrix, B∈ℛn×nu stands for the input matrix, Cr∈ℛnr×n is the output matrix, and k∈N denotes the sampling instant number. Hence, as in [7], from (1), the state predictions for the consecutive sampling instants are(2)xk+1=Axk+Buk,xk+2=A2xk+ABuk+Buk+1,⋮xk+i=Aixk+Ai−1B…ABBuk⋮uk+i−2uk+i−1,xk+i=Aixk+Ai−1B…ABB⨿1nu,N⋮⨿i−1nu,N⨿inu,Nu˜k,xk+i=Φixk+Ψiu˜k,∀i∈1,…,N,where Φi∈ℛn×n and Ψi∈ℛn×Nnu are used to represent the N-step-ahead prediction map for Linear Time Invariant (LTI) systems in a compactness form (2). Thus, system (2) can be reformulated using the following vector-matrix notation:(3)x˜k=Φxk+Ψu˜k,where(4)Φ≔Φ1Φ2⋮ΦN,Ψ≔Ψ1Ψ2⋮ΨN,and x˜k∈ℛN⋅n is the whole state trajectory of x=x1x2…xnT, u˜k∈ℛN⋅nu concatenates the computed sequence of u=u1u2…unuT, and y˜rk∈ℛN⋅nr is the output trajectory of yr=y1y2…ynrT:(5)x˜k≔xk+1xk+2⋮xk+N,u˜k≔ukuk+1⋮uk+N−1,y˜rk≔yrk+1yrk+2⋮yrk+N.Now, any candidate sequence of actionsu˜k has the corresponding future behavior of the system contained in the state trajectory x˜k and consequently the output trajectory y˜rk.Defining the projection matrix⨿in,N,(6)⨿in,N≔0n×n…0n×n︸i−1termsIn×n0n×n…0n×n︸N−iterms,where n is the vector length, N corresponds to the prediction horizon, and i stands for the desired term, then the state, control, and output vectors at a specific instant k+i can be obtained as follows:(7)xk+i=⨿in,N⋅x˜k,uk+i−1=Пinu,N⋅u˜k,yrk+i=Пin,N⋅y˜rk.


## 3. Proposed Relaxed Cost Function

### 3.1. Development

To find the best sequence of control action, in the present work, the value of the cost functionJu˜xk,w˜rk,up is defined over the prediction horizon k,k+N as follows:(8)J≔∑i=1N⨿inr,N⋅y˜r−⨿inr,N⋅w˜rQ2+⨿1nu,N⋅u˜−upR2+∑i=2N⨿inu,N−⨿i−1nu,Nu˜R2,where considering a past control action value up∈ℛnu, the last two terms are added to penalize the rate excursion of the control vector by the weighting matrix R∈ℛnu×nu. Moreover, in the first term, the pondering matrix Q∈ℛnr×nr penalizes the error between the output trajectory y˜rk∈ℛN⋅nr and the reference w˜rk∈ℛN⋅nr of the vector wr=w1w2…wnrT, which is expressed as follows:(9)w˜rk≔wrk+1wrk+2⋮wrk+N.Using notation (7) and (8),(10)J≔∑i=1Nyrk+i−wrk+iQ2+⨿1nu,N⋅u˜−upR2+∑i=2N⨿inu,N−⨿i−1nu,Nu˜R2.Moreover, substituting (1) and (2) in the first element of (10), the cost function is expressed in state-space representation:(11)J≔∑i=1NCrΨiu˜+CrΦixk−wrk+iQ2+⨿1nu,N⋅u˜−upR2+∑i=2N⨿inu,N−⨿i−1nu,Nu˜R2.Thus, to obtain the control law, the performance index (11) is rewritten in the following compact form:(12)J≔12u˜THu˜+F1xk+F2w˜r+F3upTu˜+Cst,where(13)H≔2∑i=1NΨiTCrTQCrΨi+⨿1nu,NTR⨿1nu,N+∑i=2N⨿inu,N−⨿i−1nu,NTR⨿inu,N−⨿i−1nu,N,F1≔2∑i=1NΨiTCrTQCrΦi,F2≔−2∑i=1NΨiTCrTQ⨿inr,N,F3≔−2⨿1nu,NTR,Cst≔∑i=1NCrΦixk−wrk+iQ2+upR2.Then, the control law is obtained from (12) as follows:(14)u˜optxk=−H−1F1xk+F2w˜r+F3up,and applying the first action of the best sequence control action u˜optxk, the MPC state feedback is(15)KMPCxk≔⨿1nu,N⋅u˜optxk,w˜rk,up.The MPC optimization problem (16) is solved at each sampling instant in which the measurements of the outputs and state variables are updated continuously:(16)Pxk:minu˜∈ℛNnuJu˜xk,w˜rk,up.


### 3.2. Definition of the Constraints

The systematic handling of constraints provided by predictive control strategies allows for significant improvements in performance over conventional control methodologies [30]. As in [7, 31], the sequence of the future actions u˜ cannot generally be freely chosen in ℛN⋅nu. At least, saturation constraints on the actuators have to be taken into account giving the following optimization problem Pxk:(17)Pxk:minu˜∈ℛNnuJu˜xks.t.gu˜xk≤0.The cost functionJu˜xk is minimized considering the constraints arranged in gu˜xk. Nevertheless, solving general optimization problems involving a high number of decision variables and high number of constraints is generally a very difficult task. For this reason, in the present work, based on the Inverse Distance Weighting (IDW) method [24, 32–36], the cost function (12) input and output constraints are included in Q and R matrices, where, to reduce computational complexity, Taylor series expansion is used to obtain the predicted manipulation u^p and predicted output y^p [23, 37–41]:(18)u^pk+1=uk+k−k−11!uk−uk−1k−k−1+⋯,(19)y^pk+1=yk+k−k−11!yk−yk−1k−k−1+⋯.The penalization at the outputQ∈ℛnr×nr is defined:(20)Q=q1⋯0⋮⋱000qp,where qp=eσ⋅γ,∀p∈1,…,nr, α is the tuning parameter set for a desired closed-loop performance, σ=α+y^p−wpminy^p−ypmax,y^p−ypmin, and γ=1/minu^p−upmax,u^p−upmin, which corresponds to an IDW function.The deviation between the predicted outputy^p and the reference wp is penalized by a function of distance y^p−wp. Furthermore, to avoid that the predicted process variable y^p remains close to the lower bound ypmin or the upper bound ypmax the miny^p−ypmax,y^p−ypmin term is added. Finally, γ term is included in order to increase the penalization while the predicted manipulation variable u^p is close to the upper or lower bounds, upmax and upmin, respectively.The penalization at the inputR∈ℛnu×nu is defined:(21)R=r1⋯0⋮⋱000rq,where rq=γ+δu^q−δuqp,∀q∈1,…,nu.The first term is included to avoid that the predicted manipulated variableu^p remains close to the upper or lower bounds, uqmax and uqmin, while the second term is added to reduce the large excursion between the increment of the predicted manipulated variable δu^q and the increment of the previous manipulated variable δuqp.The resulting optimization problemPxk expressed as the cost function Ju˜xk,w˜rk,up is minimized at each instant k by a sequence of future actions u˜optxk respecting all the constraints. As it is described, Q and R are defined as functionals of constraints depending on the current state xk, the desired trajectory w˜rk, and the past controlled input up. Thus, the matrices H, Fi, i=1,2,3 are computed online looking for relaxing the optimization problem.


### 3.3. Algorithm for the Proposed MPC Strategy

The algorithm used for the proposed discrete-time state-space constrained MPC strategy consists of(1)
Define the linear time invariant mathematical model of the physical system in state-space representation (equation (1))(2)
Compute the matricesΦi and Ψi,∀i∈1,…,N used to represent the N-step-ahead prediction map for LTI systems (equation (4))(3)
Estimate the manipulationu^p (equation (18)) and the output y^p (equation (19)) at the next sampling instant, using Taylor series expansion(4)
Compute the weighting matricesQ (equation (20)) and R (equation (19)) based on the proposed functions(5)
Compute the online matricesH, Fi, i=1,2,3 of the proposed performance index (equation (12))(6)
Obtain the control lawu˜optxk (equation (14)) and apply to the system the first action of the best sequence control action KMPCxk (equation (15))(7)
At the next sampling instantk≔k+1, the reference wr, the bounds ymin, ymax, umin, and umax, the measurements of the outputs yr, and the state variables xk are updated, and the MPC optimization problem is solved again (equation (17)); the algorithm goes to step 3


## 3.1. Development

To find the best sequence of control action, in the present work, the value of the cost functionJu˜xk,w˜rk,up is defined over the prediction horizon k,k+N as follows:(8)J≔∑i=1N⨿inr,N⋅y˜r−⨿inr,N⋅w˜rQ2+⨿1nu,N⋅u˜−upR2+∑i=2N⨿inu,N−⨿i−1nu,Nu˜R2,where considering a past control action value up∈ℛnu, the last two terms are added to penalize the rate excursion of the control vector by the weighting matrix R∈ℛnu×nu. Moreover, in the first term, the pondering matrix Q∈ℛnr×nr penalizes the error between the output trajectory y˜rk∈ℛN⋅nr and the reference w˜rk∈ℛN⋅nr of the vector wr=w1w2…wnrT, which is expressed as follows:(9)w˜rk≔wrk+1wrk+2⋮wrk+N.Using notation (7) and (8),(10)J≔∑i=1Nyrk+i−wrk+iQ2+⨿1nu,N⋅u˜−upR2+∑i=2N⨿inu,N−⨿i−1nu,Nu˜R2.Moreover, substituting (1) and (2) in the first element of (10), the cost function is expressed in state-space representation:(11)J≔∑i=1NCrΨiu˜+CrΦixk−wrk+iQ2+⨿1nu,N⋅u˜−upR2+∑i=2N⨿inu,N−⨿i−1nu,Nu˜R2.Thus, to obtain the control law, the performance index (11) is rewritten in the following compact form:(12)J≔12u˜THu˜+F1xk+F2w˜r+F3upTu˜+Cst,where(13)H≔2∑i=1NΨiTCrTQCrΨi+⨿1nu,NTR⨿1nu,N+∑i=2N⨿inu,N−⨿i−1nu,NTR⨿inu,N−⨿i−1nu,N,F1≔2∑i=1NΨiTCrTQCrΦi,F2≔−2∑i=1NΨiTCrTQ⨿inr,N,F3≔−2⨿1nu,NTR,Cst≔∑i=1NCrΦixk−wrk+iQ2+upR2.Then, the control law is obtained from (12) as follows:(14)u˜optxk=−H−1F1xk+F2w˜r+F3up,and applying the first action of the best sequence control action u˜optxk, the MPC state feedback is(15)KMPCxk≔⨿1nu,N⋅u˜optxk,w˜rk,up.The MPC optimization problem (16) is solved at each sampling instant in which the measurements of the outputs and state variables are updated continuously:(16)Pxk:minu˜∈ℛNnuJu˜xk,w˜rk,up.


## 3.2. Definition of the Constraints

The systematic handling of constraints provided by predictive control strategies allows for significant improvements in performance over conventional control methodologies [30]. As in [7, 31], the sequence of the future actions u˜ cannot generally be freely chosen in ℛN⋅nu. At least, saturation constraints on the actuators have to be taken into account giving the following optimization problem Pxk:(17)Pxk:minu˜∈ℛNnuJu˜xks.t.gu˜xk≤0.The cost functionJu˜xk is minimized considering the constraints arranged in gu˜xk. Nevertheless, solving general optimization problems involving a high number of decision variables and high number of constraints is generally a very difficult task. For this reason, in the present work, based on the Inverse Distance Weighting (IDW) method [24, 32–36], the cost function (12) input and output constraints are included in Q and R matrices, where, to reduce computational complexity, Taylor series expansion is used to obtain the predicted manipulation u^p and predicted output y^p [23, 37–41]:(18)u^pk+1=uk+k−k−11!uk−uk−1k−k−1+⋯,(19)y^pk+1=yk+k−k−11!yk−yk−1k−k−1+⋯.The penalization at the outputQ∈ℛnr×nr is defined:(20)Q=q1⋯0⋮⋱000qp,where qp=eσ⋅γ,∀p∈1,…,nr, α is the tuning parameter set for a desired closed-loop performance, σ=α+y^p−wpminy^p−ypmax,y^p−ypmin, and γ=1/minu^p−upmax,u^p−upmin, which corresponds to an IDW function.The deviation between the predicted outputy^p and the reference wp is penalized by a function of distance y^p−wp. Furthermore, to avoid that the predicted process variable y^p remains close to the lower bound ypmin or the upper bound ypmax the miny^p−ypmax,y^p−ypmin term is added. Finally, γ term is included in order to increase the penalization while the predicted manipulation variable u^p is close to the upper or lower bounds, upmax and upmin, respectively.The penalization at the inputR∈ℛnu×nu is defined:(21)R=r1⋯0⋮⋱000rq,where rq=γ+δu^q−δuqp,∀q∈1,…,nu.The first term is included to avoid that the predicted manipulated variableu^p remains close to the upper or lower bounds, uqmax and uqmin, while the second term is added to reduce the large excursion between the increment of the predicted manipulated variable δu^q and the increment of the previous manipulated variable δuqp.The resulting optimization problemPxk expressed as the cost function Ju˜xk,w˜rk,up is minimized at each instant k by a sequence of future actions u˜optxk respecting all the constraints. As it is described, Q and R are defined as functionals of constraints depending on the current state xk, the desired trajectory w˜rk, and the past controlled input up. Thus, the matrices H, Fi, i=1,2,3 are computed online looking for relaxing the optimization problem.


## 3.3. Algorithm for the Proposed MPC Strategy

The algorithm used for the proposed discrete-time state-space constrained MPC strategy consists of(1)
Define the linear time invariant mathematical model of the physical system in state-space representation (equation (1))(2)
Compute the matricesΦi and Ψi,∀i∈1,…,N used to represent the N-step-ahead prediction map for LTI systems (equation (4))(3)
Estimate the manipulationu^p (equation (18)) and the output y^p (equation (19)) at the next sampling instant, using Taylor series expansion(4)
Compute the weighting matricesQ (equation (20)) and R (equation (19)) based on the proposed functions(5)
Compute the online matricesH, Fi, i=1,2,3 of the proposed performance index (equation (12))(6)
Obtain the control lawu˜optxk (equation (14)) and apply to the system the first action of the best sequence control action KMPCxk (equation (15))(7)
At the next sampling instantk≔k+1, the reference wr, the bounds ymin, ymax, umin, and umax, the measurements of the outputs yr, and the state variables xk are updated, and the MPC optimization problem is solved again (equation (17)); the algorithm goes to step 3


## 4. Simulation and Results

The present work and the MPC strategies in [7, 31] are developed using the same computational platform for its evaluation. Hence, to solve the optimization problem in [7, 31], quadprog toolbox from MATLAB® is used. MATLAB® codes for the following study cases are available: MPC Matlab Files_MPE.


### 4.1. Example 1: Tape Transport System

The tape drive system consists of two reels to supply and file data. Here, the data transfer rate is proportional to the tape transport speed. Thus, the tape drive mechanism must be able to rapidly transport a fragile tape with an accurate tension regulation. Figure1 shows the schematic of a tape transport system where its components and variables involved are the tape stiffness and the damping denoted by K and D, the reel radii and the inertia represented as r and J, the motor torque constant Kt, and the viscous friction coefficient denoted by β.Figure 1
Tape transport system.Assuming there is no force loss across the head, the tape tensionT=T1=T2 [42]. Although the physical model of the process contains nonlinearities, in (22), a simplified state-space model is presented in continuous time [42–44]:(22)x˙t=Acxt+Bcut,yt=Ccxt,Ac=−Dr12J1+r22J2Dβr1J1−Kr1Kr2−Dβr2J2r1J1−βJ10−r2J20−βJ2,Bc=−DKtr1J1DKtr2J2KtJ100KtJ2,Cc=0r12r22100.The model has three states,x=Tω1ω2T, where T is the tension tape in N; meanwhile, ω1 and ω2 represent the supply and take-up reel in rad/s, respectively. Moreover, the system has two inputs u=u1u2T that represent the voltages applied to the reel motors in volts, and two outputs y=y1y2T which stand for the tape speed vrw at the read-write head in m/s and the tape tension T, respectively. The control strategy described in the present work is simulated using parameters from the tested tape system described in [42, 45, 46], whose parameters are summarized in Table 1. Considering that the motors are nominally identical, for both motors, it is used as the same motor torque constant Kt and viscous friction coefficient β [45, 46].Table 1
Tape transport system parameters.
Symbol
Parameter
Value
K
Tape stiffness
2 × 103 N/m
D
Damping
2 N s/m2
r
1
Radius of supply reel
21.2 × 10−3 m
r
2
Radius of take-up reel
9.75 × 10−3 m
J
1
Moment of inertia of the supply reel
14.2 × 10−6 Kg m2
J
2
Moment of inertia of the take-up reel
10.35 × 10−6 Kg m2
K
t
Motor torque constant
24.8 × 10−3 N m/V
β
Viscous friction coefficient
1.03 × 10−4 N m s/radDiscretizing (23) with a sampling time τ=0.1s and considering a prediction horizon N=10 with a tuning parameter α=3.7, the simulation results are shown in Figure 2. It is divided in two main parts, the first part corresponds to the outputs and the inputs of the system using the MPC with the relaxed cost function, while the second part shows the coefficient values of the weighting matrices Q and R. As it is shown, computing the coefficient values q1 and q2, the constrained tape speed y1 and tape tension y2 present a good performance. Here, y1 do not present overshoot and has a maximum settling time of 0.1 seconds while y2 has a maximum overshoot of 3% and a settling time of 0.3 seconds. On the contrary, the motor voltages u1 for the supply reel and u2 for the take-up reel remains inside their lower and upper bounds by using the computed coefficient values r1 and r2. Here, the values of r1 are lower than the values of r2, due to the rates of the manipulations.Figure 2
Simulation results for tape transport system under the MPC strategy with the relaxed cost function.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)In order to see the closed-loop stability of the system, the stability indicatorSN is defined as in [7]:(23)SN≔maxj∈1,…,nλjA−BKN,where λj stands for the eigenvalues of the system for j∈1,…,n and KN≔⨿1nu,N⋅H−1F1xk is the closed-loop gain. Figure 3 shows that the system remains stable during the test. Finally, Table 2 is presented to compare the execution time between the present work and previous works [7, 31].Figure 3
Stability behavior.Table 2
Execution time comparison with previous works.
MPC
Total (s)
Computation ofu˜k (s)
Percentage consumption (%)
[5]
2.919
1.936
66.3
[22]
2.907
2.160
74.3
This work
0.892
0.012
1.3As it can be seen, the total execution time is reduced, by taking advantage of the relaxed cost function. Here, the computation ofu˜k is 1.924 seconds and 2.148 seconds faster than the computation of the manipulations using [7, 31] MPCs strategies. Henceforth, the percentage consumption of time to obtain the control actions is notably decreased.


### 4.2. Example 2: Jet Transport Aircraft

The Jet Transport Aircraft Boeing 747 in high-lift configuration addresses complex geometries and physical phenomena that make the controller design a difficult process. Figure4 illustrates the Jet Transport Aircraft with its components and variables involved such as the angles β and φ and the angular velocities ψ and θ.Figure 4
Jet transport aircraft. (a) Top view. (b) Front view.Although the physical model of the Boeing 747 is lengthy, in (24), the simplified state-space model during cruise flight at MATCH=0.8 and H=40,000ft is presented in continuous time [47]:(24)x˙t=Acxt+Bcut,yt=Ccxt,Ac=−0.0558−0.99680.08020.04150.598−0.115−0.03180−3.050.388−0.465000.080510,Bc=0.07290−4.750.7750.15314.300,Cc=01000001.The model has four states,x=βψθφT, where β is the sideslip angle, φ stands for the bank angle, and meanwhile ψ and θ represent the yaw and roll rate, respectively. Herein, all the angles are in rad and the angular velocities in rad/s. The system has two inputs u=u1u2T: the rudder and the aileron deflections, and two outputs y=y1y2T: the yaw rate ψ and the bank angle φ.Using a sampling timeτ=0.2s, system (24) is discretized. Then, a set of changes in the output reference and a series of variations in the constraints are used to test the present work. Considering a prediction horizon N=20 and a tuning parameter α=4.5, simulation results are shown in Figure 5. Here, using the computed coefficient values q1 and q2, the yaw rate y1 has a maximum overshoot of 2.9% with a maximum settling time of 1 second, while the bank angle y2 has a maximum overshoot of 6.6% with a maximum settling time of 1.2 seconds. Moreover, with the computed coefficient values r1 and r2, the rudder deflection u1, and the aileron deflection u2 remains inside their bounds without saturation. Here, the values of r2 are greater than the values of r1, this is because the rate of the manipulation u2 is greater than the rate of the manipulation u1. On the contrary, Figure 6 shows the stability behavior during the test. As it is shown, the system remains stable.Figure 5
Simulation results for jet transport aircraft under the MPC strategy with the relaxed cost function.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)Figure 6
Stability behavior.Finally, the execution time comparison between the present work and previous works [7, 31] is shown in Table 3. As it can be seen, the total execution time is considerably reduced due to the relaxed cost function. Here, the computation of u˜k is 6.074 seconds and 3.155 seconds faster than the computation of the manipulations using [7, 31] MPCs strategies. Then, as in example 1, the percentage consumption of time to obtain the control actions is notably decreased.Table 3
Execution time comparison with previous works.
MPC
Total (s)
Computation ofu˜k (s)
Percentage consumption (%)
[5]
8.262
6.105
73.9
[22]
4.020
3.186
79.2
This work
1.808
0.031
1.7


## 4.1. Example 1: Tape Transport System

The tape drive system consists of two reels to supply and file data. Here, the data transfer rate is proportional to the tape transport speed. Thus, the tape drive mechanism must be able to rapidly transport a fragile tape with an accurate tension regulation. Figure1 shows the schematic of a tape transport system where its components and variables involved are the tape stiffness and the damping denoted by K and D, the reel radii and the inertia represented as r and J, the motor torque constant Kt, and the viscous friction coefficient denoted by β.Figure 1
Tape transport system.Assuming there is no force loss across the head, the tape tensionT=T1=T2 [42]. Although the physical model of the process contains nonlinearities, in (22), a simplified state-space model is presented in continuous time [42–44]:(22)x˙t=Acxt+Bcut,yt=Ccxt,Ac=−Dr12J1+r22J2Dβr1J1−Kr1Kr2−Dβr2J2r1J1−βJ10−r2J20−βJ2,Bc=−DKtr1J1DKtr2J2KtJ100KtJ2,Cc=0r12r22100.The model has three states,x=Tω1ω2T, where T is the tension tape in N; meanwhile, ω1 and ω2 represent the supply and take-up reel in rad/s, respectively. Moreover, the system has two inputs u=u1u2T that represent the voltages applied to the reel motors in volts, and two outputs y=y1y2T which stand for the tape speed vrw at the read-write head in m/s and the tape tension T, respectively. The control strategy described in the present work is simulated using parameters from the tested tape system described in [42, 45, 46], whose parameters are summarized in Table 1. Considering that the motors are nominally identical, for both motors, it is used as the same motor torque constant Kt and viscous friction coefficient β [45, 46].Table 1
Tape transport system parameters.
Symbol
Parameter
Value
K
Tape stiffness
2 × 103 N/m
D
Damping
2 N s/m2
r
1
Radius of supply reel
21.2 × 10−3 m
r
2
Radius of take-up reel
9.75 × 10−3 m
J
1
Moment of inertia of the supply reel
14.2 × 10−6 Kg m2
J
2
Moment of inertia of the take-up reel
10.35 × 10−6 Kg m2
K
t
Motor torque constant
24.8 × 10−3 N m/V
β
Viscous friction coefficient
1.03 × 10−4 N m s/radDiscretizing (23) with a sampling time τ=0.1s and considering a prediction horizon N=10 with a tuning parameter α=3.7, the simulation results are shown in Figure 2. It is divided in two main parts, the first part corresponds to the outputs and the inputs of the system using the MPC with the relaxed cost function, while the second part shows the coefficient values of the weighting matrices Q and R. As it is shown, computing the coefficient values q1 and q2, the constrained tape speed y1 and tape tension y2 present a good performance. Here, y1 do not present overshoot and has a maximum settling time of 0.1 seconds while y2 has a maximum overshoot of 3% and a settling time of 0.3 seconds. On the contrary, the motor voltages u1 for the supply reel and u2 for the take-up reel remains inside their lower and upper bounds by using the computed coefficient values r1 and r2. Here, the values of r1 are lower than the values of r2, due to the rates of the manipulations.Figure 2
Simulation results for tape transport system under the MPC strategy with the relaxed cost function.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)In order to see the closed-loop stability of the system, the stability indicatorSN is defined as in [7]:(23)SN≔maxj∈1,…,nλjA−BKN,where λj stands for the eigenvalues of the system for j∈1,…,n and KN≔⨿1nu,N⋅H−1F1xk is the closed-loop gain. Figure 3 shows that the system remains stable during the test. Finally, Table 2 is presented to compare the execution time between the present work and previous works [7, 31].Figure 3
Stability behavior.Table 2
Execution time comparison with previous works.
MPC
Total (s)
Computation ofu˜k (s)
Percentage consumption (%)
[5]
2.919
1.936
66.3
[22]
2.907
2.160
74.3
This work
0.892
0.012
1.3As it can be seen, the total execution time is reduced, by taking advantage of the relaxed cost function. Here, the computation ofu˜k is 1.924 seconds and 2.148 seconds faster than the computation of the manipulations using [7, 31] MPCs strategies. Henceforth, the percentage consumption of time to obtain the control actions is notably decreased.


## 4.2. Example 2: Jet Transport Aircraft

The Jet Transport Aircraft Boeing 747 in high-lift configuration addresses complex geometries and physical phenomena that make the controller design a difficult process. Figure4 illustrates the Jet Transport Aircraft with its components and variables involved such as the angles β and φ and the angular velocities ψ and θ.Figure 4
Jet transport aircraft. (a) Top view. (b) Front view.Although the physical model of the Boeing 747 is lengthy, in (24), the simplified state-space model during cruise flight at MATCH=0.8 and H=40,000ft is presented in continuous time [47]:(24)x˙t=Acxt+Bcut,yt=Ccxt,Ac=−0.0558−0.99680.08020.04150.598−0.115−0.03180−3.050.388−0.465000.080510,Bc=0.07290−4.750.7750.15314.300,Cc=01000001.The model has four states,x=βψθφT, where β is the sideslip angle, φ stands for the bank angle, and meanwhile ψ and θ represent the yaw and roll rate, respectively. Herein, all the angles are in rad and the angular velocities in rad/s. The system has two inputs u=u1u2T: the rudder and the aileron deflections, and two outputs y=y1y2T: the yaw rate ψ and the bank angle φ.Using a sampling timeτ=0.2s, system (24) is discretized. Then, a set of changes in the output reference and a series of variations in the constraints are used to test the present work. Considering a prediction horizon N=20 and a tuning parameter α=4.5, simulation results are shown in Figure 5. Here, using the computed coefficient values q1 and q2, the yaw rate y1 has a maximum overshoot of 2.9% with a maximum settling time of 1 second, while the bank angle y2 has a maximum overshoot of 6.6% with a maximum settling time of 1.2 seconds. Moreover, with the computed coefficient values r1 and r2, the rudder deflection u1, and the aileron deflection u2 remains inside their bounds without saturation. Here, the values of r2 are greater than the values of r1, this is because the rate of the manipulation u2 is greater than the rate of the manipulation u1. On the contrary, Figure 6 shows the stability behavior during the test. As it is shown, the system remains stable.Figure 5
Simulation results for jet transport aircraft under the MPC strategy with the relaxed cost function.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)Figure 6
Stability behavior.Finally, the execution time comparison between the present work and previous works [7, 31] is shown in Table 3. As it can be seen, the total execution time is considerably reduced due to the relaxed cost function. Here, the computation of u˜k is 6.074 seconds and 3.155 seconds faster than the computation of the manipulations using [7, 31] MPCs strategies. Then, as in example 1, the percentage consumption of time to obtain the control actions is notably decreased.Table 3
Execution time comparison with previous works.
MPC
Total (s)
Computation ofu˜k (s)
Percentage consumption (%)
[5]
8.262
6.105
73.9
[22]
4.020
3.186
79.2
This work
1.808
0.031
1.7


## 5. Conclusions

This paper presents a discrete-time state-space MPC approach for multivariable systems. Based on the IDW method and the concept of Taylor series expansion, a relaxed performance index with constraints defined in the online weighting matrices is proposed to compute the control action. Thus, as in study cases, the proposed MPC strategy is used to control a tape transport system and a jet transport aircraft during cruise flight.Simulation results show that the proposed MPC strategy with the relaxed cost function has a good performance, no matter abrupt changes of set-points and constraints occur, even at the same time. Additionally, compared with two well-known discrete-time state-space MPC strategies, there is a significant improvement on the execution time without affecting the tracking performance. The percentage consumption of time to compute the best sequence of control actionsu˜k is 1.3% for the tape transport system and 1.7% for the jet transport aircraft. Henceforth, it takes almost 0.1 milliseconds for the tape transport system and 0.12 milliseconds for the jet transport aircraft to obtain the manipulation uk that minimizes the proposed cost function while respecting the constraints. Thus, the proposed MPC strategy with the relaxed cost function is mainly suitable for constrained multivariable real processes with fast dynamics.



---

*Source: 7485865-2020-03-31.xml*