In this paper the solution of optimal control problems under actuator excursion and rate limits is considered. One way of addressing internal limits of a servomechanism in optimal control problems is through the introduction of input and state constraints. However, this modeling approach may introduce unnecessary conservatism in the optimal solution of the control problem. In the present work an actuator formulation based on Mixed Logical Dynamics (MLD) is developed to address this issue, eliminating the conservatism. These conditions are implemented through a Mixed-Integer Linear Programming (MILP) formulation with binary variables. The application of the formulation is illustrated using a model of a general transport aircraft. In this scenario both the standard optimal control formulation with state constraints and the novel MILP formulation proposed in this paper are used to maximize a heading change of the aircraft at a fixed terminal time. It is shown that the MILP encoding yields a lower cost function value compared to the standard optimal control formulation at the expense of greater computational resources.
? Rubens Afonso acknowledges the support of CAPES (fellowship proc. #88881.145490/2017-01) and the Federal Ministry for Education and Research of Germany through the Alexander von Humboldt Foundation. Johannes Diepolder acknowledges the support by the DFG grant HO4190/8-2. ?? Copyright c 2020 for this paper by its authors. Use permitted under Creative Commons License Attribution 4.0 International (CC BY 4.0).
Actuator saturations are ubiquitously present in control systems and limit the capacity of steering the plant. Such saturations usually originate from the physics of the operation of the actuators, e.g. the aperture of a valve, which cannot lie outside the interval between 0 and 100%. Hence, even if the controller commands an aperture of 110%, this is not physically achievable. If left unaddressed by the controller, such limitations may compromise the performance or even the stability of the closed-loop system (see e.g. [10]). Therefore, it is important to consider these limitations within the design of control laws.
In the optimal control framework, actuator rate or excursion saturations are usually modeled via reformulation of such limitations as constraints within the optimization problem, as in [6]. Therefore, the optimal control is determined so that these constraints are not violated. On the other hand, the saturations are enforced by the actuators themselves in practice, even if the controls are not determined to respect them. This effectively limits the set of candidate solutions to the optimization problem to those that do not violate the constraints a priori, whereas in reality no solution would violate these constraints by their very nature. This simplified modeling approach does not come without cost, as the limitation a priori may entail suboptimal solutions [4].
For example, consider a valve: even if the control entails an aperture of more than 100% of its capacity, such an aperture cannot be physically achieved. Furthermore, suppose that the same control command is sent to other elements in the control loop. One can hypothesize that the same command is used to drive a heater whose limits are not reached by the same command value that would entail an aperture of 110% of the valve. Therefore, by using a command that is smaller in magnitude so that the valve is not steered to open more than 100% a priori, one is also limiting the output of the heater, yielding a possibly suboptimal solution. In contrast, the saturation of the valve can be considered in a hybrid model, in which it behaves linearly between 0% and 100%, but does not leave this interval in spite of commands in that sense. Then, the command could exceed the corresponding value and the heater might be driven to its optimal state, thus enabling an optimal solution to the real control problem. In other words, if the control could be chosen freely on such constrained arcs of the more limited actuator, then it could optimally command the remaining one.
In the present paper, the optimal control problem involves deciding commands that drive actuators in closed-loop. An aircraft model is used and focus is given in the lateral motion, with the actuators being the left and right ailerons, the rudder, as well as the elevator. Each of these control surfaces is modeled as a second-order system, driven by the references of the heading as well as the normal and lateral load factors. The lateral acceleration is regulated to zero throughout the maneuver and the aim is to optimize the heading command to achieve the maximum heading change within a given fixed time interval. The excursion of the control surface actuators and their rates of change are limited. Moreover, a fault is assumed to take place in the left aileron, whereas the right one as well as the rudder and elevator operate normally. This fault results in a severely reduced excursion for the affected surface. This gives rise to the kind of problem discussed above, as the saturation is usually modeled by state constraints in the optimal control problem, leading to unnecessary conservatism. An early example of such modeling in the context of aircraft control using Linear Programming can be found in [8]. Another example can be found in [2], where the fault is a jamming of one surface, which leads to a necessity of performing a different control allocation to the remaining actuators.
In this work, the actuator saturations are modeled as physical limits that are imposed a posteriori, i.e., the excursion is saturated if the command would lead to a violation. The proposed optimization problem considers the saturation, but does not regard it as an state constraint, rather treats it using logic: whenever the limit is to be exceeded, it is saturated and the derivative is set to zero. Thus, the proposed approach is not as conservative as the usual modeling using state constraints, but yields an optimization problem that takes the saturation into account.
The presence of logical clauses in the optimization problem leads to a Mixed Logical Dynamics (MLD) model, as the control and state variables are subject to the continuous-time dynamics of the aircraft and servomechanisms, whereas the saturations introduce logical switching between linear models for the actuators. For implementing the logic, binary optimization variables are used in conjunction with the real-valued variables representing the controls and states. The presence of these binary variables entails a Mixed-Integer Linear Programming (MILP) problem, a class known to be NP-Complete [5]. Consequently, it is expected that the computational cost is increased as compared to the usual optimal control formulation. Therefore, in the proposal the number of binary variables is reduced by blocking them for a number of successive time steps, which than introduces a compromise between optimality and computational burden.
The remainder of the present paper is divided as follows. The usual formulation as an Optimal Control Problem (OCP) and the proposed MLD reformulation are presented in Section 2. The aircraft model for the illustrative example considered in this work is outlined in Section 3. Simulation results comparing the proposal and the usual OCP are shown and discussed in Section 4. Finally, concluding remarks and suggestion for future work are given in Section 5. 2 2.1
Let us first introduce an optimal control problem formulation which considers
actuator limits using purely state dependent constraints. The discrete
LinearTime-Invariant (LTI) model is defined as
xk+1 = Adxk + Bduk, k = 0, . . . , N − 1,
(
Similarly to the state constraints, the limits of all unbounded states xf,lb (xf,ub) are set to −∞ (+∞).
Considering the optimal control problem on a horizon of N time-steps the linear problem formulation is: ulb ≤ uk ≤ uub. xlb ≤ xk ≤ xub.
x0 = x(t0)
(
In the following, the equality and inequality constraints which need to be introduced for the actuator model depending on the values of the decision variables at each discrete time point are described. Let Ceq denote a set of functions which collects all equality constraints related to the actuator model. Similarly, let Ciq denote a set of functions for all inequality constraints.
Depending on the choice of the binary variables di,k, i = 1, 2, 3, 4 for each
time-step tk → tk+1 the constraints which need to be added to these sets are
described in Alg. 1. This procedure is called for all k ∈ {0, . . . , N − 1} in a
loop and updates the equality and inequality constraints in the sets Ceq and Ciq
(both initialized as empty sets before the first iteration of the loop). Note that
the anticipated rate and position states x˜p and x˜v (cf. lines 2 and 3 in Alg. 1)
are not components of the state vector. They are rather predictions without
considering the effect of the saturations. Observe that the last constraint which
is added to the set Ciq of inequality constraints
d1,k + d2,k + d3,k + d4,k − 1 ≤ 0
(
Note that all constraint functions c : Rp → R are of the general form c(zc) = acT zc − bc, with a parameter vector zc ∈ Rp, a constraint vector ac ∈ Rp as well as the scalar right-hand-side bc ∈ R. For an inequality constraint function c(zc) ≤ 0 the big-M formulation can be written as where C(d) is a condition depending on a vector of binary variables d which takes either the value true or false (respectively zero or one) and M(c(zc)) returns a
Algorithm 1 Actuator model constraints
1: procedure C(Ceq, Ciq, d1...4,k, xk, xk+1, uk, av, ap, bv, bp,xv,lb,xv,ub,xp,lb,xp,ub)
. Add actuator constraints for a time-step tk → tk+1
2: x˜v ← avTxk + bvTuk . Predicted rate state
3: x˜p ← apTxk + bpTuk . Predicted position state
4: if d1,k then
5: Ciq ← Ciq ∪ {x˜v − xv,lb ≤ 0}
6: Ceq ← Ceq ∪ {xv,k+1 − xv,lb = 0}
7: Ceq ← Ceq ∪ {xp,k+1 − apTxk − bpTuk = 0}
8: else if d2,k then
9: Ciq ← Ciq ∪ {−x˜v + xv,ub ≤ 0}
10: Ceq ← Ceq ∪ {−xv,k+1 + xv,ub = 0}
11: Ceq ← Ceq ∪ {xp,k+1 − apTxk − bpTuk = 0}
12: else if d3,k then
13: Ciq ← Ciq ∪ {x˜p − xp,lb ≤ 0}
14: Ciq ← Ciq ∪ {x˜v ≤ 0}
15: Ceq ← Ceq ∪ {xp,k+1 − xp,lb = 0}
16: Ceq ← Ceq ∪ {xv,k+1 = 0}
17: else if d4,k then
18: Ciq ← Ciq ∪ {−x˜p + xp,ub ≤ 0}
19: Ciq ← Ciq ∪ {−x˜v ≤ 0}
20: Ceq ← Ceq ∪ {−xp,k+1 + xp,ub = 0}
21: Ceq ← Ceq ∪ {xv,k+1 = 0}
22: else
23: Ceq ← Ceq ∪ {xp,k+1 − apTxk − bpTuk = 0}
24: Ceq ← Ceq ∪ {xv,k+1 − avTxk − bvTuk = 0}
25: end if
26: Ciq ← Ciq ∪ {d1,k + d2,k + d3,k + d4,k − 1 ≤ 0}
27: return Ceq,Ciq
28: end procedure
conservative big-M constant for relaxing the constraint. For this big-M constant
it needs to be ensured that given an admissible region Zc for the parameters zc
the value of the constraint function is always lower than this constant. The
minimum cost function value of the following Linear Program M(c(zc))
minimize
zc ∈ Zc
− c(zc)
yields the maximum value of the constraint function. However, it remains to
determine the admissible set Zc for the parameters zc. Observe that for all
constraints the vector zc merely contains the discretized control and state
variables. As such, the admissible set is defined by the limits of these variables.
For the controls these limits are readily available as all inputs to the system
are assumed bounded (cf. Eq. (
The extension to equality constraints is straight-forward due to the fact that an equality constraint may be expressed using two inequality constraints and the same procedure for the computation of the big-M constant can be employed for each of the inequality constraints. For completeness, the full list of constraints including the corresponding conditions which need to be introduced using the big-M notation is provided in Table 1. The dynamic model under consideration represents the closed-loop system of a generic twin-engine, fixed-wing aircraft. The plant, which is described as a rigid body, is controlled by the primary control surfaces, i.e. the left and right ailerons as well as the elevator and the rudder. The inner-loop controller for this aircraft features integral as well as proportional parts and translates a roll angle command Φc, as well as normal and lateral load factor commands, nz,c and ny,c, to desired control surface deflections. For each of the primary control surfaces (ailerons, elevator, and rudder) the dynamics are modeled by linear second-order models (natural frequency ωn = 40 rad/s, relative damping ζ¯ = 0.707) of the form d x(t) = dt x˙ (t)
0 1
−ωn −2ζ¯ωn
2
xx˙ ((tt)) + k0a δ(t),
(
It is important to mention that this particular model assumes a single actuator for the ailerons, as one considers that symmetric commands are sent to the left and right ailerons. In the present work, the left aileron is assumed to be in a fault condition, which limits its excursion severely more than the nominal one, whereas the right aileron and the rudder operate under nominal conditions. Therefore, in order to separate both ailerons, left and right deflections ξl and ξr, as well as rates ξ˙l and ξ˙r, respectively, are considered instead of a single deflection ξ and rate ξ˙. Their effect is taken into account by hypothesizing that each of these deflections generate half of the moment that would be generated by an equivalent ξ. Therefore, their coefficient is half of the coefficient of ξ and then their deflections are summed. Under nominal operation, i.e., when both ailerons are subject to the same limits, this should yield the same result as a single aileron symmetric deflection variable applied to them. 4 4.1
Longitudinal
The maneuver used to evaluate the proposal in the present paper consists of minimizing the heading Ψ of the aircraft at a fixed final instant tf = 10 s. Such a maneuver corresponds to maximizing the turn of the aircraft nose in the horizontal plane towards the clockwise sense. Constraints are imposed on the variables related to the excursions and rates of the actuators as indicated in Table 4. The commanded roll angle is also limited to avoid unrealistic reference values, thus Φc ∈ [−45, 45] deg. Moreover, the final roll angle was set to Φ(tf ) = 0 deg and the roll ωx and yaw ωz rates were subject to terminal constraints ωx(tf ) = ωz(tf ) = 0 deg/s, aiming at returning the aircraft to a leveled flight condition at the end of the maneuver.
For the discretization employed in the optimization problems the time interval [0, tf ] was sampled with N = 500 points at a regular interval. However, for the case of the proposed MILP approach, such a quantity of samples entails a large number of binary variables, which are linked to the computational demands for solving the optimization problem. On the other hand, sampling with lower frequency might compromise the capability of predicting accurately the system dynamics, especially of the fastest modes. Therefore, as a compromise solution, the binary variables were kept constant for 10 samples, enabling the real-valued variables to be sampled at a suitable frequency and at the same time reducing the number of binary decision variables to one tenth.
A time limit of 120 s was imposed for the MILP solution, i.e., if after this period the optimal solution was not yet determined, the solver returned the feasible solution with the lowest cost so far, if any was available. 4.2
The results using the standard optimal control and MILP modeling can be seen in Fig. 1 for the heading Ψ , roll Φ and sideslip β angles, as well as the commanded roll Φc angle. The achieved final values of the heading are Ψ (tf ) = −9.2 deg imposing the fault as constraints and Ψ (tf ) = −28.4 deg with the MILP approach. This difference can be attributed to a less conservative usage of the roll angle with the usual modeling approach, since it does not achieve values lower than −20 deg, whereas with the MILP approach the minimum value is lower than −40 deg. The sideslip angle is kept low throughout the maneuver with both approaches.
Figure 2 depicts the actuator deflections and rates. Both approaches respect the constraints. However, the MILP scheme enables an asymmetric usage between the left (faulty) and right (healthy) ailerons. Indeed, the right aileron is driven to its lower limit of −20 deg between 0 and 2 s. In contrast, without the MILP formulation, the right aileron is steered only between −2 and 2 deg. It is important to remark that the MILP approach preserved the constraints imposed on the left aileron, remaining between −2 and 2 deg. Moreover, one can see that the availability of more actuator authority in the right aileron leads the MILP solution to use the rudder more intensively, as larger sideslip angle variations are seen with the more aggressive curve flown. 4.3
For validation purposes, the inputs obtained from the MILP solution were applied to the nonlinear model around the trimmed condition. It is important to remark that this model, besides being nonlinear, involves also the complete dynamics of the aircraft, including the longitudinal movement. The results are presented and discussed in the present subsection.
The heading achieved in the nonlinear simulation can be seen in Fig. 3 to be lower in absolute value in the nonlinear simulation, reaching Ψ (tf ) = −26.1 deg, a reduction of 2.3 deg in the span as compared to the linear simulation. However, this remains substantially higher in absolute value when compared to the one obtained when the MILP formulation was not used. Regarding the roll angle, no considerable difference can be noticed in Fig. 3. On the other hand, the sideslip angle can be seen to reach slightly higher absolute values in the nonlinear simulation, albeit still acceptable.
In Figure 4, it can be seen that the left aileron deflection remains at the lower bound for a longer period in the nonlinear simulation as compared to the results with the linear model. Moreover, the right aileron deflection follows the same profile obtained in the linear simulation, with small differences due to different amplitudes and durations of the oscillations in the rate associated to this aileron. 0.5 -0.5 -1 0 1 1 1 2 2 2 3 3 3 4 4 4 5 5 5 6 6 6 7 7 7 8 8 8 9
10 9
10 9 10 Saturations as constraints
MILP model Fig. 1. Commands and outputs treating saturation as constraints and with the proposed approach. 8
10 8 10 20 10 0 -10 -20 20 0 -20 20 0 0 0 -20 0 2 2 2 2 0 -2
0 20 0 -20
0 10 0 -10 0 2 2 2 4 4 4 6 6 6 8 10 8 10 Saturations as constraints
MILP model Fig. 2. Actuator deflections and rates treating saturation as constraints and with the proposed approach. 4 4 4 6 6 6 8
10 8 10
The same applies for the rudder. Overall, one can ascertain that the nonlinear simulation confirms the applicability of the controls determined by the proposal and their superior performance in terms of the minimization of the heading.
0 -10 -20 -30 20
0 -20 -40 1 0 -1 -2 0 0 0 1 1 1 2 2 2 3 3 3 4 4 4 5 5 5 6 6 6 7 7 7 8 8 8 9
10 9
10 9 10 MILP linear simulation
MILP nonlinear simulation Fig. 3. Commands and outputs with the proposed approach simulated with the linear and nonlinear models.
2 2 2 4 4 4 6 6 6 8 8 8 10 10 10 20 0 -20 20 0 -20 20 0 0 0 -20
0 MILP linear simulation
MILP nonlinear simulation 2 0 -2
0 20 0 -20 0 5 0 -5 -10 0 2 2 2 4 4 4 6 6 6 8 8 8 10 10 10
In the present work a MLD formulation for modeling internal limits (rate and position) for actuation systems was developed. Under this approach the logical conditions were implemented through a MILP formulation with binary variables. The additional degree of freedom enabled by this formulation results in an improved solution compared to the standard optimal control formulation in the simulation example of the lateral movement of a general transport aircraft. On the other hand, the proposed MILP encoding yields a higher computational burden.
Future work could involve increasing the computational efficiency for solving the resulting MILP problem e.g. through the techniques proposed in [9, 11] for reducing the number of binary variables in the MILP encoding.