# Offline Model Predictive Control Based on Weighted Projection over Polytopes

Authors: Fabio Andrade de Almeida
Journal: Journal of Applied Mathematics
 (2015)

---

## Abstract

This work presents a novel offline model predictive control technique for tracking of constrained systems. The quadratic programming problem, commonly found in constrained control methods, is replaced by sequential offline set projections based on priority given to the decision variables. If a preference is established in terms of which decision variables are more desirable, the optimization problem can be solved by sequentially choosing the most important variables and performing a membership test with the projection of the constraint closed set over the related dimensions. Thus, real-time optimization is replaced by offline projection operations and online one-dimensional membership tests. This concept of decision variable prioritization is then applied to a form of model predictive control: feasible target tracking. Three quadratic programming problems are replaced by the proposed method. In the first problem, attainable steady-state demands are computed based on the performance of the plant. The reachable target command is then filtered in terms of dynamic admissibility, creating feasible inputs to the plant. Finally, the control is computed considering the current state and disturbance vectors along with the feasible and attainable command. Simulations of the method executing a path-following task are presented, demonstrating its benefits with negligible online computational burden.


---

## Body

## 1. Introduction

Model predictive control has been a matter of extensive research in recent decades because it provides a straightforward way to address the control of constrained multivariate systems [1]. The computed control to be applied in the plant is such that it satisfies the system constraints; therefore, an online constrained optimal control must be computed at every time instant. Though being initially conceived for applications in the process industry, other engineering fields are taking advantage of the developed techniques, such as in space systems [2] and aeronautical guidance and control [3, 4]. The facility of MPC to deal with constraints and its common appearance with traditional linear-quadratic controllers makes this technique suitable for the design of fault-tolerant flight control systems and also to allow operations near vehicle operational bounds.Despite the fact that successful flight tests have already been executed to perform constrained control of fast dynamical modes of an aircraft [5], the related computational load is still high for the present generation of flight control computers. However, in the last decade, new approaches have been presented to decrease the computational burden of the online control calculations. Explicit MPC [6] has become the main strategy for moving offline most of the constrained control law computations. In the other direction, new methods were recently proposed [7, 8] to improve numerical performance during real-time execution.Nevertheless, the online computation of the control actions depends mostly on linear or quadratic programming solvers, which do not offer guarantees of finite-time execution because of their iterative procedures. This makes MPC a noncertifiable technique for civil aviation control software, where deterministic algorithms are mandatory [9]. Explicit MPC could be an alternative for this problem, but it is still restricted for control of low-order systems and for regulation purposes [10]. An attempt to perform robust variable tracking, instead of nominal regulation, can increase the complexity of the state-space formulation in a way that makes it hard or even impossible to compute offline solutions.Therefore, the main contribution of this work is the proposal of a novel method to move offline the nondeterministic computations of MPC, leaving only set membership verifications during real-time operation. This is achieved by proposing a new form to solve quadratic convex optimization problems and applying it on the feasible target tracking model predictive control (FTT-MPC) method [11], producing a new form of offline MPC. Based on the assumption of the existence of priority among controlled outputs and manipulated variables (inputs), the traditional quadratic programming problem is replaced by offline computation and projection of polyhedral sets. This new model predictive control method is then applied to perform airspeed, flight path angle, and turn rate tracking of a small UAV model during a path-following task with constraints on angle of attack, throttle, and aerodynamic surfaces deflections as well as with considerable plant-model mismatches, demonstrating adequate flight envelope estimation and protection.The structure of this paper is organized as follows. Section2 establishes the constrained control problem and initial assumptions. Section 3 presents a new method to compute the projection of a point over a polytope, based on prioritization of the decision vector components. Section 4 is the main part of this work, where the feasible target tracking MPC is converted to an offline implementation. Finally, Section 5 applies the technique to the longitudinal and lateral control of an unmanned aircraft, showing a practical and potential application of the proposed MPC method.


## 2. Problem Formulation

Let the discrete-time state-space model be defined by(1)xk+1=Φxk+Γuk+Γddk,dk+1=dk,yk=Exk,zk=Hyk,
where xk∈Rn is the state vector, uk∈Rm is the control vector, yk∈Rp is the vector of observed variables, dk∈Rq is the constant vector of disturbances, and zk∈Rq, q≤m, is the vector of controlled variables. The disturbance vector is included to model mismatches between the plant and the nominal model as well as constant external disturbances acting on the plant. Note that the number of disturbance variables is equal to the number of controlled outputs q.The objective of the control system is to asymptotically eliminate the tracking error given a piecewise constant reference commandrc, which is
(2)zk→k→∞rc
in the presence of a disturbance dk and given constraints on the state and control vectors xk∈X,uk∈U, where X and U are closed, bounded, and convex sets expressed by linear inequalities. In case rc is not achievable, a reachable target rss will be computed with some priority among the controlled variables. Additionally, the dynamical protection characteristic of the controller maintains the outputs inside their admissible ranges by computing feasible target values of zss.Before proceeding with the assumptions, some useful definitions are given. A polyhedron is the intersection of a finite number of closed halfspaces:(3)P=v∈Rn∣Lcv≤kc,
where Lc and kc are, respectively, a matrix and a vector which represent linearly the set of halfspaces. A polytope is a bounded polyhedron. Given a polyhedron P⊂V×Y, where V and Y are subspaces, that is, sets with lower dimensions that P, the projection π of P over V, is defined as
(4)πVP=v∈V∣∃y∈Y,v,y∈P,
and thesliceσ of P at a vector r∈Rp is defined as
(5)σrP=v∈Rn∣v,r∈P,Assumption 1.
The pairΦ,E is assumed to be detectable with E full row rank. The following condition is equivalent:
(6)rankΦ-IΓE0=n+p.Assumption 2.
It is assumed that the state and disturbance vectors are estimated by a nominally asymptotically stable observer. Becausedk∈Rq, the observer proposed by Maeder et al. [12] is chosen for this work. Other choices of observers [13] would increase the dimension of the disturbance vector and add additional numerical complexity to the controller synthesis.Assumption 3.
The disturbance is bounded and lies in a compact convex polyhedron containing the origin in its interior.


## 3. Weighted Projection of a Point on a Polytope

In the next section, it will be presented a new form of MPC that basically shapes input references and computes the constrained control, by projecting points on convex polyhedral sets. The projected point is commonly computed by solving a quadratic programming (QP) problem [14] and several solvers are available to calculate the solution; however, they comprise iterative procedures with no guarantee of finite time execution. This section presents a novel and simple way to project a point over a polytope without need of online optimization.The Euclidean orthogonal projection ofv0 on a polytope is given by the solution of the QP [14]:
(7)min⁡vJv0=v-v0Tv-v0subjectto:v∈V.
We also refer to oblique 2-norm projection of the minimizer of the following weighted cost function:
(8)min⁡yJy0=y-y0TWy-y0subjectto:y∈Y,
where W∈Rn×n is a positive definite diagonal matrix. Clearly, both formulations are equivalent by doing the variable transformation vi=wiyi, i=1,…,n. Here the weights are chosen to influence the solution, giving relative different emphasis to some components in the minimization of the distance between the optimal solution y* and y0.It is proposed to solve (8) with a special choice of each diagonal element of the W matrix: the largest individual weight penalizes the most important component in terms of distance minimization, followed by the second component, and so on. Giving enough scale separation between the diagonal elements, it is possible to solve (8) through sequential projections and slices of polytopes.The proposed prioritized computation ofyp* is given as follows. Let us consider the feasible set Y. Its projection on the most important component y1 creates the one-dimensional set Y1. Thus, the prioritized value y1p* is obtained by verifying whether y10∈Y1 or not. If it is, y1p*=y10. If not, y1p* will be equal to the nearest extremum of Y1. Next, the set Y is projected over the two most important components, generating Y2. Because y1p* has been already selected, the choice of y2p* is done over the slice of Y2 at y1p*, in which it is also a one-dimensional set. This sequence of projections and slices continues until the last component yn.Figure1 exemplifies and compares, using a two-dimensional case, two solutions of (8), one considering W=I, and the other with weighting prioritization, that is, w1≫w2. The established setpoint (y10,y20) is not achievable because it lies outside the admissible set. In the unweighted projection computation, the pair (y1*,y2*) is the closest point to (y10,y20) and is perpendicular to the related active constraint. However, considering y1 as the preferred variable to attain, the demand prioritization technique computes a solution where y2p* is the result of the prioritized choice y1p*. Its selection is simply performed by a membership test with the projection of the two-dimensional polytope on the first dimension, as shown in the lower part of Figure 1.Figure 1
Orthogonal and prioritized projections of a point on a polytope.The extension of the proposed method to compute decision variables in higher-dimensional spaces and its application to the problem of constrained control must take into consideration the definition of the projection operation. As stated by (4), the projection of a high-dimensional decision space on the most important variables subspace guarantees the existence of the remaining decision elements. Doing sequential projections clearly decreases the sets of the variables with low priority, even with degeneration to a point, but they will not be empty.The presented method of solving (8) is particularly helpful for online applications. It is worth observing that the successive projections of the set Y can be done completely offline. It only remains to calculate online the sequence of slices, which are computationally inexpensive.


## 4. Offline Feasible Model Predictive Control

While input constraints are relatively straightforward to manipulate, output constraints must be mapped onto the input space through some sort of dynamical prediction. Thus, the input applied to the system will consider the future evolution of the constrained outputs. Even basic formulations of MPC naturally compute optimal inputs that force the outputs to stay inside their limitations. However, not all reference values to be tracked are numerically feasible; that is, the imposed limitations and the system demand may not be compatible. In addition, even if feasibility is found, closed-loop stability may not be assured because the constraints are active not only at the time instant of the input computation, but also at all future instants.Previous work presented the feasible target tracking MPC technique [11], a three-stage control method where each part comprises a quadratic programming problem. In the present work, the offline feasible model predictive control method (OFF-MPC) is proposed, replacing each optimization of the FTT-MPC by prioritized projections on polytopes that comprises the relevant system constraints.The offline feasible MPC technique, schematically shown in Figure2, guarantees both feasibility and closed-loop stability by first computing statically reachable commands, which respects the complete set of constraints when the system reaches the steady-state condition. The disturbance vector of (1) is required to correct the steady-state output limitations, that is, the admissible envelope, when differences between the model and the real system occur. After this first computation, dynamically feasible commands are calculated, to assure that the input and output constraints will be satisfied during the transient. In this step, both state and disturbance vectors are necessary because the dynamically feasible demand depends on the state at each time instant. Finally, the optimal constrained control is obtained in the MPC, considering a feasible and reachable demand.Figure 2
Overall scheme of the offline feasible MPC.


### 4.1. Sets of Statically Admissible Commands

The separation between the determination of steady-state target values and constrained control computation creates distinct optimization problems, which are individually solvable through the demand prioritization technique. Thus, in a first step, the polyhedral set of input, state, and output constraints that act in the steady-state is created and projected offline over the highest priority variable to track, creating an interval of admissible and reachable values. In the sequence, another projection is performed over the second most important variable, and this operation repeats until the last variable, whose computation is simply performed by online slicing the complete constraint set with the choice of prioritized variables.The equality constraints represented by the state-space and target equations can be converted into inequality constraints through introduction of arbitrarily small tolerance vectorsϵx and ϵz. Hence, (1) can be stated, at steady state, as
(9)Φ-Ixss+Γuss+Γddk≤ϵxI,Φ-Ixss+Γuss+Γddk≥-ϵxI,HExss-rss≤ϵzI,HExss-rss≥-ϵzI,uss∈λU,xss∈λX,
which defines the polyhedron Pss spanned by the augmented vector xssTxssTrssTdkTT. The projection of Pss over rssTdkTT generates a polyhedron Pw that encodes statically admissible commands rss to the real plant, provided that the disturbance is known or estimated. The scalar λ∈(0,1) slightly shrinks the constraints at steady state to allow the set computations of the next subsection.If there is no priority among references to be tracked, the statically admissible vectorrss is the orthogonal projection of rc on the slice of Pw at estimated values of dk. However, the prioritization and ordering of each element of rc give the possibility of avoiding online quadratic programming computation, being replaced by an offline sequence of polyhedra projections.The first required offline operation is a projection ofPw over r1ssdkTT to create the polyhedron Pw1. The procedure follows with the projection of Pw over r1ssr2ssdkTT to compute Pw2. This is performed by all ordered elements of rc.The online slice ofPw1 with estimated values of dk gives the admissible range of r1ss. Thus, in the sequence, slices are performed with the disturbance and selected attainable reference values. This sequence of offline projections and online slices is followed until the last demand rq, as shown by Algorithm 1.Algorithm 1:Sets of statically admissible commands.
Off-line computation
(1) Sort the reference vectorrc in order of importance
(2) Create the polyhedronPss
(3) ProjectPss over [rssTdkT]T to create the polyhedron Pw
(4) Fori = 1 to q
(a) ProjectPw over [r1ss⋯rissdkT]T to create Pwi
On-line computation
(1) SlicePw1 at estimated values of dk
(2) Check ifr1 lies inside the admissible range of r1ss
(3) Fori = 2 to q
(a) SlicePwi at selected [r1ss⋯r(i-1)ss]T and estimated dk
(b) Check ifri lies inside the admissible range of riss


### 4.2. Sets of Dynamically Admissible Commands

Once the attainable targetrss is computed, the next step of the OFF-MPC method is the computation of the dynamically admissible target zss, which guarantees the existence of a feasible control action to steer states and controls to the target values xss∈λX and uss∈λU. The starting point is the proposed control law predicted for N steps [15]:
(10)uj=-Kxj-xss+uss+cjj=k,…,k+N-1,uj=-Kxj-xss+ussj≥k+N.
The model predictive controller calculates the contribution cj over N steps. A reparameterization that relates the target state vector with the target control vector and disturbance is also proposed:
(11)xss=I-Φ-1Γuss+I-Φ-1Γddk.
Therefore, assuming that rss is piece-wise constant, it will be convenient to extend the state-space vector as xa=xjTdjTussTrssTT, leading to the following extended linear model:
(12)xj+1dj+1ussrss=Φ-ΓKΘdΓdΓΘ00I0000I0000Ixjdjussrss+Γ000cj,
where Θ=K(I-Φ)-1Γ+I and Θd=ΓK(I-Φ)-1+I. In a similar way of that established in the previous section, the relation between the statically admissible demand rss, steady-state control uss, and the disturbance dk=dss is found through the formulation of the following constraint:
(13)HEI-Φ-1Γddk+HEI-Φ-1Γuss-rss≤ϵrI-HEI-Φ-1Γddk-HEI-Φ-1Γuss+rss≤-ϵrI
in addition to state limitations xj∈X. Because the reference value rss is statically admissible, uss fully respects the control limitations and does not need to be accounted for again. For stability purposes, it is assumed that, after N steps, the extended state vector xa at j=k+N lies inside an invariant set O∞. This extension of the state-space vector is similar to those proposed by Chisci and Zappa [16] and to recent techniques of economic MPC [17].One possible way of interpreting the role of the MPC control is that the inputscjj=kk+N-1 are responsible for steering the extended state in N that moves into the invariant set O∞. The set of all states that can be steered to O∞ after N steps is the domain of attraction O∞(N) of the complete closed-loop constrained system. This set is computed by first performing recursive substitutions of the extended dynamic system into the constraints, creating the set O∞c and subsequently projecting it over xak to eliminate the MPC contributions in the domain of attraction formulation.Here the priority among the demands again eliminates the need for online optimization. Performing a projection ofO∞c on the subspace spanned by xjTdjTr1ssT removes the dependency on uss, creating the polyhedron Pr1. The slice of Pr1 with estimated values of xk and dk gives the dynamically admissible range of r1ss. It is worth noting the similarity with the procedure adopted to obtain the static set in the previous subsection, but instead of considering only the estimated disturbance vector, the current state vector is also considered here to guarantee transient constraint satisfaction.Once the most important demand is filtered in a dynamic sense and compared to the admissible range, a new commandz1ss is created. The sequence of dynamic filtering, which is summarized in Algorithm 2, follows the same steps of the static admissibility.Algorithm 2:Sets of dynamically admissible commands.
Off-line computation
(1) Create invariant setO∞ of the extended
closed-loop state-space
(2) Create the setO∞c by recursion
(3) Fori = 1 to q
(a) ProjectO∞c over [xjTdjTr1ss⋯riss]T to create Pri
On-line computation
(1) SlicePr1 at estimated values of xk and dk
to create the one-dimensional setZ1ss
(2) Check ifr1ss is inside Z1ss to obtain z1ss
(3) Fori = 2 to q
(a) SlicePri at selected [z1ss⋯z(i-1)ss]T and
known values ofxk and dk to create
the one-dimensional setZiss
(b) Check ifriss is inside Ziss to obtain zissRemark 4.
The reparameterization presented before is suitable for both square and nonsquare systems, whereq≤m and I-Φ admits an inverse. If one or both conditions do not hold, the steady-state reparameterization suggested by Shead and Rossiter [18] is applicable, and the extension of the state space is performed considering the null space of Φ-IΓ, instead of uss.


### 4.3. Control Calculation

The computation ofzss introduces one additional degree of freedom in the constrained controller. Because the dynamically filtered demand is calculated from O∞c, the state and control vectors at steady state respect the related constraints. The feedforward control input uss is calculated as
(14)uss=HEI-Φ-1Γ-1zss-HEI-Φ-1Γddk.Assumption 5.
There are no elements of the state vector integrating elements of the controlled vector. The following condition is equivalent:(15)rankΦ-IΓHE0=n+q.Assumption 6.
(1) LetΨ=ΨT∈Rm×m be a diagonal positive definite matrix.
(2) LetK∈Rm×n be a feedback gain such that Φ-ΓK is Hurwitz.
The constrained control problem is to determine the optimal control sequenceCN*=cjj=kk+N-1 that minimizes, at each time instant k, the following cost function:
(16)Jk=∑j=kk+N-1‍cjTΨcj
subject to cj∈slice(O∞c,[x^kTd^kTzssT]T), j∈(k,k+N-1).Theorem 7.
Assume that Assumptions1 through 6 hold. Assume that the piecewise constant command rss is statically admissible. Thus, for any initial estimated state, disturbance and command (x^k,d^k,rss)∈O∞(N), the dynamically admissible command zss and control law uk=κ(x^k,d^k,rss) steer the system to the statically admissible setpoint rss, while respecting the constraints.Proof.
First, the control lawuk=-K(x^k-xss)+uss+ck is a feasible solution for the constrained control problem because xak∈O∞(N). Let us define the optimal cost function at the sample instant k as Jk*=∑j=kk+N-1‍cjTΨcj with feasible control sequence CN*=cjj=kk+N-1. At the following sample time k+1, define the control sequence CN-1=cjj=k+1k+N-1 and also Jk+1=∑j=k+1k+N‍cjTΨcj. Thus, Jk*-Jk+1=ckTΨck. Considering that at k+1 the optimal cost is Jk+1*, it follows that Jk+1≥Jk+1* and therefore Jk+1*-Jk*≤ckTΨck.
Hence, the optimal cost sequenceJk*k=0∞ can never increase. Moreover, if Jk*=Jk+1*, then ckTΨck=0, which means that the extended state-space xak already lies already in the invariant set O∞, from where the unconstrained linear feedback law uj+N=-K(xj+N-xss)+uss steers the controlled vector to rss when k→∞.Remark 8.
The fact thatrss is piecewise constant means that the above stability result holds, only if changes on rss occur after the system reaches its equilibrium around steady state. In this case, stability can be assured, but not when the reference is modified while the system dynamics is still prevailing.The calculation of the MPC contributions can also follow the concept of decision variable prioritization, ordering the elementsc1,c2,…,cm. From the polyhedron O∞c, one can project it offline over variables xjTdjTrssTc1T to create the set Pc1. The known values of xk,dk, and rss=zss are then used to online slice this set and provide an admissible range for c1. The procedure is repeated to generate the sets Pc1,…,Pcm and is summarized in Algorithm 3. Note that the MPC contribution will be zero unless it is required to steer the states to the invariant set O∞. Additionally, note that it is necessary to compute only the first m elements of CN* because the remaining inputs are not applied in the plant at time instant k. Thus, the online computational burden does not depend on the control horizon size, which affects only the offline projections of O∞c on Algorithms 2 and 3.Algorithm 3:Control calculation.
Off-line computation
(1) Sort the control contribution vectorck
(2) Fori = 1 to m
(a) ProjectO∞c over [xjTdjTrssTc1⋯cm]T to create Pci
On-line computation
(1) Calculatexss and uss from the dynamically
admissible vectorzss
(2) SlicePc1 at estimated values of xk, dk and rss=zss
(3) Check ifc1=0 lies inside the admissible values
(4) Fori = 2 to m
(a) SlicePci at computed [c1⋯c(i-1)]T and at
values ofxk, dk and rss=zss
(b) Check ifci=0 lies inside the admissible values
(5) Computeuk=-K(xk-xss)+uss+ck


### 4.4. Optimality of the Prioritized Constrained Solution

During normal operation, inside the system limitations, the control consists of a regulation around state and control steady-state values. This regulation is governed by the feedback linear gainK. Thus, optimality is easily achieved if K is obtained through optimal control methods (e.g., LQR).Given a statically admissible reference to track, the controller offers the computation of bothzss and ck to respect system constraints. First, if some constraint becomes active, the dynamical filter is the first subsystem to act (Algorithm 2), computing zss in such a way that the MPC is capable of steering the states inside the extended invariant set after N moves. The ordered choice of the elements of zss does not affect the size of the domain of attraction of the controller because there are MPC nonlinear contributions ck to maintain the system stable while tracking. This is assured by the projection operation on O∞c.The ordering of controlled variables seems suitable for multivariable systems where there are interdependence among demands. This type of relationship can occur when operating at the limits of the system performance, what it is exactly MPC does, or following a desired optimal trajectory computed by some kind of higher level optimization layer. For example, the optimal climb schedule for an aircraft generally relates flight path angle (the angle that the airspeed vector makes with the ground) with the airspeed. Thus, rate or angle of climb is first selected, followed by the selection of airspeed. Depending on the aircraft performance, the airspeed is unique for a given altitude change demand. The example of Section5 will explore this aspect, presenting the prioritization of demands as a mean to achieve better aircraft climb performance.After the first filter, the contributionck is calculated (Algorithm 3) considering the previous filtered values of xss and uss. The existence of nonempty sets Pci, i=1 to m, is guaranteed by the computation of zss. Again, the prioritization is only a method to select the contributions inside each set Pci. Therefore, the way that each constrained decision variable is computed (orthogonal or prioritized projection) does not enlarge or reduce either the sets Ziss, i=1 to q, or Pci.Any choice of the MPC correction insidesliceO∞c,x^kTd^kTzssTT is such that keeps the system inside the limitations. The minimization of the cost function given by (16) provides an optimal choice among all possible solutions. If the cost function prioritizes some variables among others, through selection of individual diagonal elements of Ψ, the presented technique solves the problem of oblique projection of a point on a set, finding an optimal solution for the problem stated in (16). Based on what is depicted in Figure 1, it is possible that the prioritized projection leads, in the control inputs with less priority, to larger ck components than those computed through orthogonal unweighted projection. Thus, the ordering of control inputs makes more use of the last elements to maintain the closed-loop system inside its limitations, since the first components of CN* have preference to be set closer or equal to zero.


## 4.1. Sets of Statically Admissible Commands

The separation between the determination of steady-state target values and constrained control computation creates distinct optimization problems, which are individually solvable through the demand prioritization technique. Thus, in a first step, the polyhedral set of input, state, and output constraints that act in the steady-state is created and projected offline over the highest priority variable to track, creating an interval of admissible and reachable values. In the sequence, another projection is performed over the second most important variable, and this operation repeats until the last variable, whose computation is simply performed by online slicing the complete constraint set with the choice of prioritized variables.The equality constraints represented by the state-space and target equations can be converted into inequality constraints through introduction of arbitrarily small tolerance vectorsϵx and ϵz. Hence, (1) can be stated, at steady state, as
(9)Φ-Ixss+Γuss+Γddk≤ϵxI,Φ-Ixss+Γuss+Γddk≥-ϵxI,HExss-rss≤ϵzI,HExss-rss≥-ϵzI,uss∈λU,xss∈λX,
which defines the polyhedron Pss spanned by the augmented vector xssTxssTrssTdkTT. The projection of Pss over rssTdkTT generates a polyhedron Pw that encodes statically admissible commands rss to the real plant, provided that the disturbance is known or estimated. The scalar λ∈(0,1) slightly shrinks the constraints at steady state to allow the set computations of the next subsection.If there is no priority among references to be tracked, the statically admissible vectorrss is the orthogonal projection of rc on the slice of Pw at estimated values of dk. However, the prioritization and ordering of each element of rc give the possibility of avoiding online quadratic programming computation, being replaced by an offline sequence of polyhedra projections.The first required offline operation is a projection ofPw over r1ssdkTT to create the polyhedron Pw1. The procedure follows with the projection of Pw over r1ssr2ssdkTT to compute Pw2. This is performed by all ordered elements of rc.The online slice ofPw1 with estimated values of dk gives the admissible range of r1ss. Thus, in the sequence, slices are performed with the disturbance and selected attainable reference values. This sequence of offline projections and online slices is followed until the last demand rq, as shown by Algorithm 1.Algorithm 1:Sets of statically admissible commands.
Off-line computation
(1) Sort the reference vectorrc in order of importance
(2) Create the polyhedronPss
(3) ProjectPss over [rssTdkT]T to create the polyhedron Pw
(4) Fori = 1 to q
(a) ProjectPw over [r1ss⋯rissdkT]T to create Pwi
On-line computation
(1) SlicePw1 at estimated values of dk
(2) Check ifr1 lies inside the admissible range of r1ss
(3) Fori = 2 to q
(a) SlicePwi at selected [r1ss⋯r(i-1)ss]T and estimated dk
(b) Check ifri lies inside the admissible range of riss


## 4.2. Sets of Dynamically Admissible Commands

Once the attainable targetrss is computed, the next step of the OFF-MPC method is the computation of the dynamically admissible target zss, which guarantees the existence of a feasible control action to steer states and controls to the target values xss∈λX and uss∈λU. The starting point is the proposed control law predicted for N steps [15]:
(10)uj=-Kxj-xss+uss+cjj=k,…,k+N-1,uj=-Kxj-xss+ussj≥k+N.
The model predictive controller calculates the contribution cj over N steps. A reparameterization that relates the target state vector with the target control vector and disturbance is also proposed:
(11)xss=I-Φ-1Γuss+I-Φ-1Γddk.
Therefore, assuming that rss is piece-wise constant, it will be convenient to extend the state-space vector as xa=xjTdjTussTrssTT, leading to the following extended linear model:
(12)xj+1dj+1ussrss=Φ-ΓKΘdΓdΓΘ00I0000I0000Ixjdjussrss+Γ000cj,
where Θ=K(I-Φ)-1Γ+I and Θd=ΓK(I-Φ)-1+I. In a similar way of that established in the previous section, the relation between the statically admissible demand rss, steady-state control uss, and the disturbance dk=dss is found through the formulation of the following constraint:
(13)HEI-Φ-1Γddk+HEI-Φ-1Γuss-rss≤ϵrI-HEI-Φ-1Γddk-HEI-Φ-1Γuss+rss≤-ϵrI
in addition to state limitations xj∈X. Because the reference value rss is statically admissible, uss fully respects the control limitations and does not need to be accounted for again. For stability purposes, it is assumed that, after N steps, the extended state vector xa at j=k+N lies inside an invariant set O∞. This extension of the state-space vector is similar to those proposed by Chisci and Zappa [16] and to recent techniques of economic MPC [17].One possible way of interpreting the role of the MPC control is that the inputscjj=kk+N-1 are responsible for steering the extended state in N that moves into the invariant set O∞. The set of all states that can be steered to O∞ after N steps is the domain of attraction O∞(N) of the complete closed-loop constrained system. This set is computed by first performing recursive substitutions of the extended dynamic system into the constraints, creating the set O∞c and subsequently projecting it over xak to eliminate the MPC contributions in the domain of attraction formulation.Here the priority among the demands again eliminates the need for online optimization. Performing a projection ofO∞c on the subspace spanned by xjTdjTr1ssT removes the dependency on uss, creating the polyhedron Pr1. The slice of Pr1 with estimated values of xk and dk gives the dynamically admissible range of r1ss. It is worth noting the similarity with the procedure adopted to obtain the static set in the previous subsection, but instead of considering only the estimated disturbance vector, the current state vector is also considered here to guarantee transient constraint satisfaction.Once the most important demand is filtered in a dynamic sense and compared to the admissible range, a new commandz1ss is created. The sequence of dynamic filtering, which is summarized in Algorithm 2, follows the same steps of the static admissibility.Algorithm 2:Sets of dynamically admissible commands.
Off-line computation
(1) Create invariant setO∞ of the extended
closed-loop state-space
(2) Create the setO∞c by recursion
(3) Fori = 1 to q
(a) ProjectO∞c over [xjTdjTr1ss⋯riss]T to create Pri
On-line computation
(1) SlicePr1 at estimated values of xk and dk
to create the one-dimensional setZ1ss
(2) Check ifr1ss is inside Z1ss to obtain z1ss
(3) Fori = 2 to q
(a) SlicePri at selected [z1ss⋯z(i-1)ss]T and
known values ofxk and dk to create
the one-dimensional setZiss
(b) Check ifriss is inside Ziss to obtain zissRemark 4.
The reparameterization presented before is suitable for both square and nonsquare systems, whereq≤m and I-Φ admits an inverse. If one or both conditions do not hold, the steady-state reparameterization suggested by Shead and Rossiter [18] is applicable, and the extension of the state space is performed considering the null space of Φ-IΓ, instead of uss.


## 4.3. Control Calculation

The computation ofzss introduces one additional degree of freedom in the constrained controller. Because the dynamically filtered demand is calculated from O∞c, the state and control vectors at steady state respect the related constraints. The feedforward control input uss is calculated as
(14)uss=HEI-Φ-1Γ-1zss-HEI-Φ-1Γddk.Assumption 5.
There are no elements of the state vector integrating elements of the controlled vector. The following condition is equivalent:(15)rankΦ-IΓHE0=n+q.Assumption 6.
(1) LetΨ=ΨT∈Rm×m be a diagonal positive definite matrix.
(2) LetK∈Rm×n be a feedback gain such that Φ-ΓK is Hurwitz.
The constrained control problem is to determine the optimal control sequenceCN*=cjj=kk+N-1 that minimizes, at each time instant k, the following cost function:
(16)Jk=∑j=kk+N-1‍cjTΨcj
subject to cj∈slice(O∞c,[x^kTd^kTzssT]T), j∈(k,k+N-1).Theorem 7.
Assume that Assumptions1 through 6 hold. Assume that the piecewise constant command rss is statically admissible. Thus, for any initial estimated state, disturbance and command (x^k,d^k,rss)∈O∞(N), the dynamically admissible command zss and control law uk=κ(x^k,d^k,rss) steer the system to the statically admissible setpoint rss, while respecting the constraints.Proof.
First, the control lawuk=-K(x^k-xss)+uss+ck is a feasible solution for the constrained control problem because xak∈O∞(N). Let us define the optimal cost function at the sample instant k as Jk*=∑j=kk+N-1‍cjTΨcj with feasible control sequence CN*=cjj=kk+N-1. At the following sample time k+1, define the control sequence CN-1=cjj=k+1k+N-1 and also Jk+1=∑j=k+1k+N‍cjTΨcj. Thus, Jk*-Jk+1=ckTΨck. Considering that at k+1 the optimal cost is Jk+1*, it follows that Jk+1≥Jk+1* and therefore Jk+1*-Jk*≤ckTΨck.
Hence, the optimal cost sequenceJk*k=0∞ can never increase. Moreover, if Jk*=Jk+1*, then ckTΨck=0, which means that the extended state-space xak already lies already in the invariant set O∞, from where the unconstrained linear feedback law uj+N=-K(xj+N-xss)+uss steers the controlled vector to rss when k→∞.Remark 8.
The fact thatrss is piecewise constant means that the above stability result holds, only if changes on rss occur after the system reaches its equilibrium around steady state. In this case, stability can be assured, but not when the reference is modified while the system dynamics is still prevailing.The calculation of the MPC contributions can also follow the concept of decision variable prioritization, ordering the elementsc1,c2,…,cm. From the polyhedron O∞c, one can project it offline over variables xjTdjTrssTc1T to create the set Pc1. The known values of xk,dk, and rss=zss are then used to online slice this set and provide an admissible range for c1. The procedure is repeated to generate the sets Pc1,…,Pcm and is summarized in Algorithm 3. Note that the MPC contribution will be zero unless it is required to steer the states to the invariant set O∞. Additionally, note that it is necessary to compute only the first m elements of CN* because the remaining inputs are not applied in the plant at time instant k. Thus, the online computational burden does not depend on the control horizon size, which affects only the offline projections of O∞c on Algorithms 2 and 3.Algorithm 3:Control calculation.
Off-line computation
(1) Sort the control contribution vectorck
(2) Fori = 1 to m
(a) ProjectO∞c over [xjTdjTrssTc1⋯cm]T to create Pci
On-line computation
(1) Calculatexss and uss from the dynamically
admissible vectorzss
(2) SlicePc1 at estimated values of xk, dk and rss=zss
(3) Check ifc1=0 lies inside the admissible values
(4) Fori = 2 to m
(a) SlicePci at computed [c1⋯c(i-1)]T and at
values ofxk, dk and rss=zss
(b) Check ifci=0 lies inside the admissible values
(5) Computeuk=-K(xk-xss)+uss+ck


## 4.4. Optimality of the Prioritized Constrained Solution

During normal operation, inside the system limitations, the control consists of a regulation around state and control steady-state values. This regulation is governed by the feedback linear gainK. Thus, optimality is easily achieved if K is obtained through optimal control methods (e.g., LQR).Given a statically admissible reference to track, the controller offers the computation of bothzss and ck to respect system constraints. First, if some constraint becomes active, the dynamical filter is the first subsystem to act (Algorithm 2), computing zss in such a way that the MPC is capable of steering the states inside the extended invariant set after N moves. The ordered choice of the elements of zss does not affect the size of the domain of attraction of the controller because there are MPC nonlinear contributions ck to maintain the system stable while tracking. This is assured by the projection operation on O∞c.The ordering of controlled variables seems suitable for multivariable systems where there are interdependence among demands. This type of relationship can occur when operating at the limits of the system performance, what it is exactly MPC does, or following a desired optimal trajectory computed by some kind of higher level optimization layer. For example, the optimal climb schedule for an aircraft generally relates flight path angle (the angle that the airspeed vector makes with the ground) with the airspeed. Thus, rate or angle of climb is first selected, followed by the selection of airspeed. Depending on the aircraft performance, the airspeed is unique for a given altitude change demand. The example of Section5 will explore this aspect, presenting the prioritization of demands as a mean to achieve better aircraft climb performance.After the first filter, the contributionck is calculated (Algorithm 3) considering the previous filtered values of xss and uss. The existence of nonempty sets Pci, i=1 to m, is guaranteed by the computation of zss. Again, the prioritization is only a method to select the contributions inside each set Pci. Therefore, the way that each constrained decision variable is computed (orthogonal or prioritized projection) does not enlarge or reduce either the sets Ziss, i=1 to q, or Pci.Any choice of the MPC correction insidesliceO∞c,x^kTd^kTzssTT is such that keeps the system inside the limitations. The minimization of the cost function given by (16) provides an optimal choice among all possible solutions. If the cost function prioritizes some variables among others, through selection of individual diagonal elements of Ψ, the presented technique solves the problem of oblique projection of a point on a set, finding an optimal solution for the problem stated in (16). Based on what is depicted in Figure 1, it is possible that the prioritized projection leads, in the control inputs with less priority, to larger ck components than those computed through orthogonal unweighted projection. Thus, the ordering of control inputs makes more use of the last elements to maintain the closed-loop system inside its limitations, since the first components of CN* have preference to be set closer or equal to zero.


## 5. Simulation Example

The test aircraft used for the simulations is a six-degree-of-freedom model of the unmanned aerial vehicle Acauã, which also takes into account actuator and engine dynamics [19]. Acauã is a propeller-driven fixed-wing unmanned testbed used at the Brazilian Aeronautics and Space Institute (IAE). Linear models for control synthesis have been obtained with airspeed of 27 m/s and a pressure altitude of 620 m, at zero sideslip.The simulation scenario consists of following a path defined by geodetic waypoints and related courses, as shown by Table1. In this work, two-dimensional Dubins paths [20, 21] were used to connect each waypoint. Briefly, a 2D Dubins path is a path that flies over a certain waypoint on its desired course makes a right or left turn followed by a straight line and finishes with a right or left turn, flying over the next waypoint on a certain course. The selected turn radius for the circular transitions was 190 m, which corresponds to a nominal rate of turn of 8.14 deg/s at an airspeed of 27 m/s. The on-off turn rate commands created by 2D Dubins paths can practically be considered as piecewise constant references to be followed by the MPC-based autopilot, which is favorable in terms of closed-loop stability (Theorem 7).Table 1
List of waypoints.
Number
Lat. (deg.)
Long. (deg.)
Course (deg.)
1
-
21.9842
-
47.3426
0
2
-
21.9765
-
47.3474
120
3
-
21.9765
-
47.3377
240Thus, the lateral guidance system [22] will be responsible for creating autopilot demands to steer the vehicle straight or through right or left turns with a constant radius. In the vertical channel, piecewise constant flight path angle demands were generated to explore the full flight envelope of the aircraft model. Figure 3 shows schematically the complete closed-loop system for the simulation.Figure 3
Closed-loop system for UAV trajectory control.


### 5.1. Longitudinal Constrained Autopilot

Consider the unmanned aerial vehicle continuous-time linearized longitudinal model described in [19]:
(17)Along=-0.18490.9732-0.1182-9.8067-0.016-2.4980.967300.0227-33.97-5.924-00010,Blong=3.582-0.7128-0.0008-0.1970.0012-41.1700.
The state vector is composed of the true airspeed (m/s), angle of attack (rad), pitch rate (rad/s), and the pitch angle (rad), respectively. The control vector is composed of the engine throttle and the elevator angular deflection. Full state observation (E=I) is assumed, and it is intended to control both airspeed and flight path angle without offset with priority given to the second variable. Thus,
(18)Hlong=10000-101.
The feedback gain K was obtained by first defining an LQR continuous-time cost function and subsequently performing the discretization [23]. The elements of the Q matrix in the LQR cost function were chosen to penalize deviations of airspeed and flight path angle; thus, Q=HlongTdiag13300Hlong. The R matrix was set to R=diag100820. The dynamical system and the weight matrices were discretized with Ts= 0.05 s, from which the matrices Φ and Γ are obtained, along with the feedback gain K. The disturbance model was chosen to have two variables, and the related distribution matrix was set to Γd=Γ. The observer gains were computed using the linear Kalman method.The inputs are constrained as follows:δtc∞≤0.5 and δec∞≤10 deg. The disturbances are assumed to respect the bounds d1∞≤0.5, d2∞≤0.05. It is required to protect the aircraft against aerodynamic stall while maneuvering, imposing α∈(-8.08,3.36) deg, which corresponds to a total admissible interval of (-3.06,8.38) deg. The constraints at steady state were scaled through λ=0.99.The setsPw1, Pw2 (Algorithm 1); Pr1, Pr2 (Algorithm 2); and Pc1 and Pc2 (Algorithm 3) were computed offline through the appropriate numerical routine [24] in the MATLAB environment, setting ϵx=ϵz=ϵr=10-4.2 and the control horizon N=2.


### 5.2.  Lateral Constrained Autopilot

The same values of the control horizonN, tolerances ϵz,ϵr, and the sampling time from the longitudinal controller were used in the synthesis of the constrained lateral autopilot. The state-space continuous-time linear lateral model is given by
(19)Alat=-0.21230.0868-0.98680.3618-5.6541-1.73110.465406.8358-0.0697-0.5304001.00000.08750,Blat=00.08586.13821.1299-0.2630-5.043600,
where the state vector is composed of the sideslip (rad), roll rate (rad/s), yaw rate (rad/s), and the bank angle (rad). Aileron and rudder angular deflections (rad) define the control vector. The track priority is given first to maintain zero sideslip and subsequently to follow the turn rate. Hence,
(20)Hlat=00101000.
After discretization of the state-space formulation, the discrete lateral gain is given by
(21)Klat=0.93020.2528-0.10690.56710.38710.0666-1.88540.1874.
Full state feedback (E=I) is considered. Constraints on control angular displacements were defined as δa∞≤10 deg and δr∞≤10 deg. In addition, the related sets of Algorithms 2 and 3 were computed offline, with disturbance vector elements constrained to d1∞≤0.05 and d2∞≤0.05.It is important to mention that Assumptions1 through 6 are fully respected in the two-controller synthesis. Assumptions 1 and 2 are accomplished by the choice of the estimator (Kalman) and full state observation. Bounds for the disturbance vector are established (Assumption 3). Finally, the choice of the controlled vector in both autopilots, that is, no states as pure integrators of controlled variables and the LQR feedback gains, guarantees Assumptions 5 and 6.


### 5.3.  Simulation Analysis

Two cases were considered for the simulation. First, nominal engine power and longitudinal stability were considered in the computation of the controllers. The second case, however, simulates an underpowered Acauã (60% of maximum nominal power) with 50% of the original static stability, with the nominal controllers. Figure4 shows the horizontal path followed by the aircraft in both cases. During the flyover of waypoints 2 and 3, an abrupt modification in the turn direction was required, but the vehicle strictly followed the computed Dubins path due to the high-gain lateral controller. The maximum lateral deviation during the two cases was 8.1 meters but with immediate guidance and control regulation to zero. The feasible commands zss and aircraft responses in the first case are shown in Figure 5. On two occasions, aggressive roll reversals were demanded by the Dubins path and guidance from the lateral autopilot, but fast and stable tracking of the turn rate was achieved with small variations of the sideslip angle.Figure 4
Horizontal path view.Figure 5
Controlled outputs and state variables in the nominal simulation.In addition to the lateral tracking task, an aggressive variation of the flight path angle was simulated by an arbitrarily large command. Regardless of the airspeed command, the highest attainable climb angle (higher priority) was selected. The angle of attack has changed between maximum and minimum values but never exceeded the limitations. The full exploration of the flight envelope for climbing and descent is shown in Figure6. Engine throttle is demanded at its maximum and minimum values, and the elevator deflection reached its minimum deflection twice. An adequate lateral path-following performance was achieved through a high-gain lateral feedback control, which manipulated aileron and rudder deflection to their limits as shown also in Figure 6. Thus, unconstrained control strategies would most likely fail to track the considered turn rates.Figure 6
Inputs in the nominal simulation.The ability to cope with plant-model mismatches is demonstrated in the case of an underpowered vehicle with reduced static stability. Figure7 shows related feasible demands and aircraft responses given the same path to be followed and flight path angle command. With only 60% of the nominal maximum power, the aircraft is not able to achieve the previous nominal climb angles, but due to the slice of Pw1 and Pw2 at the estimated disturbance values, real attainable commands are online computed. It is worth noting the demanded reduction of airspeed, which is crucial in the task of tracking the maximum climb angle. The longitudinal controller maintained performance and stability despite the fact that the real aircraft was less stable. Figure 8 shows that the elevator deflection required to maneuver was noticeably less than in the first case. Additionally, the reduced available power is shown in the evolution of the engine throttle. Again, the key aspect of offset-free tracking is the disturbance vector of the longitudinal autopilot, as shown in Figure 9. While the disturbance was practically zero, that is, a small plant-model mismatch, in the first simulation, the mismatch is present but properly taken into account in this case.Figure 7
Controlled outputs and state variables in the simulation with reduced power and stability.Figure 8
Inputs in the simulation with reduced power and stability.Figure 9
Estimated disturbances of the longitudinal controller in the simulation with reduced power and stability.


### 5.4. Aspects of Computational Burden and Optimality

The elimination of online iterative numerical routines provided by OFF-MPC reflects the drastic runtime reduction. Using the MATLAB commandtic toc, the complete constrained control computation (online sections of Algorithms 1, 2, and 3) required, on average, 347 μs in both vertical and lateral autopilots. The simulations were performed on a PC with a 3.20 GHz Intel Core i5 CPU and 4 GB of RAM. A previous work [11] presented larger runtimes, in the range of milliseconds. Recently, Hartley et al. [2] presented numerical simulations of MPC with execution times in the range of microseconds, but using iterative solvers that demand specific initialization tuning. Although the offline set projections are based on iterative numerical routines, the online segment is only a direct sequence of algebraic operations, completely eliminating the dependency on iteration stop criteria. It is worth mentioning that this numerical performance was achieved using only MATLAB-coded routines, without previous compilation that would speed up the constrained computation.Another comparison is done with the purpose to demonstrate the optimality achieved when the problem stated by (8) is solved with offline projections and online slices. Figure 10 shows engine throttle (left) and elevator deflection (right) MPC corrections during a capture of 10 deg of flight path angle, using the proposed offline method and an active-set quadratic programming solver [25] to obtain the optimal solution of (16). Note the large nonlinear correction required to maintain the elevator deflection inside the limits.Figure 10
Nonlinear MPC corrections during capture of flight path angle.Although the weighting matrixΨ was set equal to identity, both numerical methods computed identical solutions. This means that the proposed technique of computing the projection of a point on a polytope is capable of providing the same results as those calculated via quadratic programming solvers. Also, by design, elevator is the aerodynamic control surface to provide fast vertical change in flight path. This type of analytical priority is captured by both methods, with advantage to the OFF-MPC, which offers the optimal solution using less computational time. In this simulation of 6 seconds, the active-set QP required, on average, 3.5 ms, against 283 μs demanded by OFF-MPC.


## 5.1. Longitudinal Constrained Autopilot

Consider the unmanned aerial vehicle continuous-time linearized longitudinal model described in [19]:
(17)Along=-0.18490.9732-0.1182-9.8067-0.016-2.4980.967300.0227-33.97-5.924-00010,Blong=3.582-0.7128-0.0008-0.1970.0012-41.1700.
The state vector is composed of the true airspeed (m/s), angle of attack (rad), pitch rate (rad/s), and the pitch angle (rad), respectively. The control vector is composed of the engine throttle and the elevator angular deflection. Full state observation (E=I) is assumed, and it is intended to control both airspeed and flight path angle without offset with priority given to the second variable. Thus,
(18)Hlong=10000-101.
The feedback gain K was obtained by first defining an LQR continuous-time cost function and subsequently performing the discretization [23]. The elements of the Q matrix in the LQR cost function were chosen to penalize deviations of airspeed and flight path angle; thus, Q=HlongTdiag13300Hlong. The R matrix was set to R=diag100820. The dynamical system and the weight matrices were discretized with Ts= 0.05 s, from which the matrices Φ and Γ are obtained, along with the feedback gain K. The disturbance model was chosen to have two variables, and the related distribution matrix was set to Γd=Γ. The observer gains were computed using the linear Kalman method.The inputs are constrained as follows:δtc∞≤0.5 and δec∞≤10 deg. The disturbances are assumed to respect the bounds d1∞≤0.5, d2∞≤0.05. It is required to protect the aircraft against aerodynamic stall while maneuvering, imposing α∈(-8.08,3.36) deg, which corresponds to a total admissible interval of (-3.06,8.38) deg. The constraints at steady state were scaled through λ=0.99.The setsPw1, Pw2 (Algorithm 1); Pr1, Pr2 (Algorithm 2); and Pc1 and Pc2 (Algorithm 3) were computed offline through the appropriate numerical routine [24] in the MATLAB environment, setting ϵx=ϵz=ϵr=10-4.2 and the control horizon N=2.


## 5.2.  Lateral Constrained Autopilot

The same values of the control horizonN, tolerances ϵz,ϵr, and the sampling time from the longitudinal controller were used in the synthesis of the constrained lateral autopilot. The state-space continuous-time linear lateral model is given by
(19)Alat=-0.21230.0868-0.98680.3618-5.6541-1.73110.465406.8358-0.0697-0.5304001.00000.08750,Blat=00.08586.13821.1299-0.2630-5.043600,
where the state vector is composed of the sideslip (rad), roll rate (rad/s), yaw rate (rad/s), and the bank angle (rad). Aileron and rudder angular deflections (rad) define the control vector. The track priority is given first to maintain zero sideslip and subsequently to follow the turn rate. Hence,
(20)Hlat=00101000.
After discretization of the state-space formulation, the discrete lateral gain is given by
(21)Klat=0.93020.2528-0.10690.56710.38710.0666-1.88540.1874.
Full state feedback (E=I) is considered. Constraints on control angular displacements were defined as δa∞≤10 deg and δr∞≤10 deg. In addition, the related sets of Algorithms 2 and 3 were computed offline, with disturbance vector elements constrained to d1∞≤0.05 and d2∞≤0.05.It is important to mention that Assumptions1 through 6 are fully respected in the two-controller synthesis. Assumptions 1 and 2 are accomplished by the choice of the estimator (Kalman) and full state observation. Bounds for the disturbance vector are established (Assumption 3). Finally, the choice of the controlled vector in both autopilots, that is, no states as pure integrators of controlled variables and the LQR feedback gains, guarantees Assumptions 5 and 6.


## 5.3.  Simulation Analysis

Two cases were considered for the simulation. First, nominal engine power and longitudinal stability were considered in the computation of the controllers. The second case, however, simulates an underpowered Acauã (60% of maximum nominal power) with 50% of the original static stability, with the nominal controllers. Figure4 shows the horizontal path followed by the aircraft in both cases. During the flyover of waypoints 2 and 3, an abrupt modification in the turn direction was required, but the vehicle strictly followed the computed Dubins path due to the high-gain lateral controller. The maximum lateral deviation during the two cases was 8.1 meters but with immediate guidance and control regulation to zero. The feasible commands zss and aircraft responses in the first case are shown in Figure 5. On two occasions, aggressive roll reversals were demanded by the Dubins path and guidance from the lateral autopilot, but fast and stable tracking of the turn rate was achieved with small variations of the sideslip angle.Figure 4
Horizontal path view.Figure 5
Controlled outputs and state variables in the nominal simulation.In addition to the lateral tracking task, an aggressive variation of the flight path angle was simulated by an arbitrarily large command. Regardless of the airspeed command, the highest attainable climb angle (higher priority) was selected. The angle of attack has changed between maximum and minimum values but never exceeded the limitations. The full exploration of the flight envelope for climbing and descent is shown in Figure6. Engine throttle is demanded at its maximum and minimum values, and the elevator deflection reached its minimum deflection twice. An adequate lateral path-following performance was achieved through a high-gain lateral feedback control, which manipulated aileron and rudder deflection to their limits as shown also in Figure 6. Thus, unconstrained control strategies would most likely fail to track the considered turn rates.Figure 6
Inputs in the nominal simulation.The ability to cope with plant-model mismatches is demonstrated in the case of an underpowered vehicle with reduced static stability. Figure7 shows related feasible demands and aircraft responses given the same path to be followed and flight path angle command. With only 60% of the nominal maximum power, the aircraft is not able to achieve the previous nominal climb angles, but due to the slice of Pw1 and Pw2 at the estimated disturbance values, real attainable commands are online computed. It is worth noting the demanded reduction of airspeed, which is crucial in the task of tracking the maximum climb angle. The longitudinal controller maintained performance and stability despite the fact that the real aircraft was less stable. Figure 8 shows that the elevator deflection required to maneuver was noticeably less than in the first case. Additionally, the reduced available power is shown in the evolution of the engine throttle. Again, the key aspect of offset-free tracking is the disturbance vector of the longitudinal autopilot, as shown in Figure 9. While the disturbance was practically zero, that is, a small plant-model mismatch, in the first simulation, the mismatch is present but properly taken into account in this case.Figure 7
Controlled outputs and state variables in the simulation with reduced power and stability.Figure 8
Inputs in the simulation with reduced power and stability.Figure 9
Estimated disturbances of the longitudinal controller in the simulation with reduced power and stability.


## 5.4. Aspects of Computational Burden and Optimality

The elimination of online iterative numerical routines provided by OFF-MPC reflects the drastic runtime reduction. Using the MATLAB commandtic toc, the complete constrained control computation (online sections of Algorithms 1, 2, and 3) required, on average, 347 μs in both vertical and lateral autopilots. The simulations were performed on a PC with a 3.20 GHz Intel Core i5 CPU and 4 GB of RAM. A previous work [11] presented larger runtimes, in the range of milliseconds. Recently, Hartley et al. [2] presented numerical simulations of MPC with execution times in the range of microseconds, but using iterative solvers that demand specific initialization tuning. Although the offline set projections are based on iterative numerical routines, the online segment is only a direct sequence of algebraic operations, completely eliminating the dependency on iteration stop criteria. It is worth mentioning that this numerical performance was achieved using only MATLAB-coded routines, without previous compilation that would speed up the constrained computation.Another comparison is done with the purpose to demonstrate the optimality achieved when the problem stated by (8) is solved with offline projections and online slices. Figure 10 shows engine throttle (left) and elevator deflection (right) MPC corrections during a capture of 10 deg of flight path angle, using the proposed offline method and an active-set quadratic programming solver [25] to obtain the optimal solution of (16). Note the large nonlinear correction required to maintain the elevator deflection inside the limits.Figure 10
Nonlinear MPC corrections during capture of flight path angle.Although the weighting matrixΨ was set equal to identity, both numerical methods computed identical solutions. This means that the proposed technique of computing the projection of a point on a polytope is capable of providing the same results as those calculated via quadratic programming solvers. Also, by design, elevator is the aerodynamic control surface to provide fast vertical change in flight path. This type of analytical priority is captured by both methods, with advantage to the OFF-MPC, which offers the optimal solution using less computational time. In this simulation of 6 seconds, the active-set QP required, on average, 3.5 ms, against 283 μs demanded by OFF-MPC.


## 6. Conclusion

The protection of the closed-loop system against excursions outside relevant limitations is an important task that must be performed by modern control systems. Moreover, performance limitations must also be taken into account to avoid commands with unattainable outputs. In this paper, a technique was presented to perform tracking on constrained systems. A new offline model predictive control technique was proposed based on the prioritization of quadratic programming decision variables. In this way, a sequence of constraint set projections over incremental subspaces of the decision variables creates polyhedral sets that, after substitution of the previously computed decision variables, degenerates to a one-dimensional set, where the choice of the optimal decision variable turns into a simple membership verification.The concept of decision variable prioritization was applied to convert the feasible target tracking model predictive control technique into an offline constrained control method. First, whether the reference vector to be tracked is statically admissible is determined. If the demand is beyond the performance limitations of the plant, an attainable solution is calculated starting with the highest priority reference, and the remaining are determined through substitutions in the incremental projections of the static constraint set.In the sequence, the statically admissible demand is filtered, taking into consideration the feasibility of the MPC problem. Extending the state space with the reference to be followed and the disturbance model, offline sequential projections of the domain of attraction of the controller and online substitutions generate one-dimensional feasible sets from which dynamically feasible target outputs are easily computed. Finally, the nonlinear MPC contributions are calculated following the same method.The complete offline MPC method was demonstrated to perform tracking of trajectory variables of an unmanned aerial vehicle model, with constraints on throttle, aerodynamic surfaces, and angle of attack. The flight envelope was estimated and protected even when plant-model mismatches were considered, yielding offset-free tracking and a complete satisfaction of constraints.



---

*Source: 741348-2015-02-18.xml*