# An Efficient Algorithm for Solving Model Predictive Control of Switched Affine Systems

S. Leirens†, J. Buisson†, P. Bastard‡ and J. L. Coullon*

†Supélec-IETR, Hybrid System Control Group

Avenue de la Boulaie, BP 81127, F-35511 Cesson-Sévigné cedex, France

E-mail: sylvain.leirens@supelec.fr

‡Supélec, Power Systems Department

3 rue Joliot-Curie, F-91192 Gif-sur-Yvette cedex, France

*Areva T&D

9 rue Ampère, F-91345 Massy, France

**Abstract** - Switched Affine Systems (SAS) involve control inputs of mixed nature. Thus Model Predictive Control (MPC) of SAS requires suitable algorithms to compute the optimal control inputs which are applied to the system at each time step. An efficient algorithm is developed in this paper which consists in finding a first suboptimum and applying a branch cutting strategy to prune the tree of switching possibilities. As a result, the size and the number of Quadratic Programming (QP) subproblems to solve is considerably reduced. The proposed algorithm is illustrated on voltage stability control of a small power system example.

**Keywords**— Predictive control, Mixed optimization, Switched affine systems, Power systems.

## I. INTRODUCTION

Most of the literature about advanced control involves optimality aspects. Thus a lot of problem formulations include something of kind: 'find the control which minimizes the cost function...' and call for suitable optimization methods. To deal with complex constrained control problems, Model Predictive Control (MPC) has become a standard in the process industry although the on-line optimization required can become computationally expensive. To deal with hybrid systems which involve variables of mixed nature, the optimization methods needed by the MPC algorithm are confronted to combinatorial aspects. Using the Mixed Logical Dynamical (MLD) framework, the system to control is modeled by discrete time dynamical equations and a set of inequality constraints involving real and binary variables [1]. The application of MPC to hybrid systems when using this framework requires the solving of a Mixed Integer Linear or Quadratic Program (MILP/MIQP) problem which is very time consuming. This optimization method involves Branch and Bound (B&B) and Linear Programming (LP) or Quadratic Programming (QP) algorithms. For Piecewise Affine (PWA) systems whose equivalence with MLD has been proved under mild assumptions in [2], some research has been engaged in the way of reducing the combinatorial complexity. In [3], the concepts of

reachable set, one step controllable set and State Transition Graph (STG) are used. In this kind of approach, the state and input space is divided into polyhedral partitions where the system evolution is described by an affine model. Switched Affine Systems (SAS) are a subclass of PWA systems that involve either discrete or continuous valued control inputs. From the modeling point of view, a change in discrete valued control inputs can be viewed as a switching between models. The switching possibilities which can be represented as a tree shape lead to strong combinatorial problems for control.

In this paper, an efficient algorithm (which does not belong to B&B kind of algorithm) is proposed to address the mixed optimization problem induced by the application of MPC to SAS. The basic ideas of this algorithm are: first to find a first suboptimum and secondly to prune the tree of switching possibilities by cutting the branches that are not able to lead to the optimum. Thus the proposed approach leads to reduce considerably the number and the size of QP subproblems to solve.

Current research involves developing models, theory, and computational techniques for hybrid systems, see [4]. Nevertheless the use of a hybrid approach to deal with the power system voltage stability issue is quite recent. In the framework of voltage stability control, power systems can be modeled using SAS. However power systems are strongly non-linear and non-linear MPC method is limited to those problems for which a globally optimal solution can be found [5]. Previous research has provided a suitable off-equilibrium linearizing procedure which allows the application of MPC to voltage stability control of power systems [6].

The paper is organized as follows: in section II, SAS and MPC concepts are introduced and the mixed optimization issue is addressed. In section III the Mixed Pruning algorithm is presented and developed. An application example coming from the area of voltage stability control of power systems is addressed in IV and conclusions are given in section V.
## II. PROBLEM FORMULATION

The class of systems which are addressed in this paper are defined by the following equations:

$$ 
\begin{aligned}
x(k+1) &= A_i x(k) + B_i u_c(k) + f_i \\
y(k) &= C_i x(k) + D_i u_c(k) + g_i
\end{aligned} 
\quad (1) 
$$

and are commonly called discrete time SAS.

These systems own two kinds of control inputs, either continuous or discrete valued. Continuous valued ones are denoted as $u_c$ whereas discrete valued ones are denoted as $u_d$ and the couple $(u_c, u_d)$ is the complete input vector which is denoted as $u$. A particular combination of the discrete valued control inputs defines an operating mode $i$ in which the system stands.

Let $x(k) \in \mathbb{R}^n$, $u_c(k) \in \mathbb{R}^m$ and $y(k) \in \mathbb{R}^p$ denote the state, continuous valued input and output vector, respectively. Matrices of (1) have appropriate dimensions: $A_i \in \mathbb{R}^{n \times n}$, $B_i \in \mathbb{R}^{n \times m}$, $C_i \in \mathbb{R}^{p \times n}$, $D_i \in \mathbb{R}^{p \times m}$, $f_i \in \mathbb{R}^n$ and $g_i \in \mathbb{R}^p$ with $i \in \mathbb{N}$. The pair $(A_i, B_i)$ is assumed to be controllable.

### A. Model Predictive Control

MPC consists in a control algorithm based on solving an on-line optimization problem using a receding horizon approach. In [7], this is summarized in the following steps:

* At time $k$ and for the current state $x(k)$ an open-loop optimal control problem over a future time interval is on-line solved taking into account the current and future constraints.

* The first step in the optimal control sequence so obtained is applied and then the procedure is repeated at time $(k+1)$ using the current state $x(k+1)$.

The solution is converted into a closed-loop strategy by using the measured value of $x(k)$ (or else an estimate provided by some form of observer) as the current state. Stability of the resultant feedback system can be established by using the fact that the cost function acts as a Lyapunov function for the closed-loop system [8].

Consider the following model:

$$ x(l+1) = f(x(l), u(l)) \quad (2) $$

and the cost function:

$$ J_N(x(k), u_N) = \sum_{l=k}^{k+N-1} L(x(l), u(l)) \quad (3) $$

where $L(\cdot) \in \mathbb{R}^+$ and

$$ u_N = \{u(k), u(k+1), \dots, u(k+N-1)\} \quad (4) $$

is the control sequence over the prediction horizon of length $N$.

In the particular setting of affine systems and quadratic cost functions, the MPC problem can be formulated by considering the model (1) and the following cost function:

$$ J_N(x(k), u_N) = \sum_{l=k}^{k+N-1} \|y(l) - y^{ref}(l)\|_{Q_y} + \|u(l)\|_{Q_u} \quad (5) $$

where $y^{ref}$ is the output reference trajectory. 2-Norm is used, i.e. $\|v\|_Q = v^T Q v$. $Q_y$ and $Q_u$ are weighing matrices such that $Q_y \ge 0$ and $Q_u > 0$.

The MPC algorithm needs to solve the optimization problem:

$$ J_N^o(x(k)) = \min_{u_N} J_N(x(k), u_N) \quad (6) $$

subject to model constraints (1) as well as current and future input and state constraints:

$$ 
\begin{aligned}
u(l) &\in \mathbb{U} \\
x(l) &\in \mathbb{X}
\end{aligned} 
\quad (7) 
$$

where $\mathbb{U}$ and $\mathbb{X}$ are the set of constraints corresponding to $u$ and $x$.

### B. Mixed Optimization

A particular feature of the optimization problem previously formulated is that it has a mixed nature. The presence of either continuous ($u_c$) or discrete ($u_d$) valued control inputs leads to consider the searched optimal control sequence as a sequence of continuous ($u_{cN}$) as well as discrete ($u_{dN}$) valued control inputs over the prediction horizon. Then the optimization problem (6) can be written as follows:

$$ J_N^o(x(k)) = \min_{u_{dN}} \left( \min_{u_{cN}} J_N(x(k), (u_{cN}, u_{dN})) \right). \quad (8) $$

Given a discrete valued control sequence $u_{dN}$,

$$ J_N^*(x(k), u_{dN}) = \min_{u_{cN}} J_N(x(k), (u_{cN}, u_{dN})) \quad (9) $$

is the optimal cost found by solving the corresponding QP subproblem. The associated continuous valued control sequence denoted as $u_{cN}^*$ is optimal regarding to $u_{dN}$.

To find the optimal control sequence over the prediction horizon, all the possible discrete valued control
<figure><img src="image_1.png" /><figcaption>Fig. 1. Brute force strategy for mixed optimization</figcaption></figure>

sequences have to be enumerated. Figure 1 shows a tree-shaped graphical representation of this brute force strategy. The depth of the tree is fixed by the length of the receding horizon whereas the width of the tree is fixed by the number of possibilities for the discrete valued control inputs (branching factor). For each discrete valued control sequence, i.e. at each leaf of the tree (represented by a small circle on the figure), a QP subproblem still remains to be solved in order to find the associated optimal continuous control sequence. Then the sequence of discrete and continuous valued control inputs which minimizes the cost function (5) is the searched optimal control sequence.

The advantage of this strategy is that it guarantees to find the searched optimum but unfortunately, and it is its main drawback, it leads to a combinatorial explosion with a possible huge number of QP subproblems to solve. Moreover it should be noted that QP subproblem solving is time consuming, of course depending on the size of this subproblem.

### III. MIXED PRUNING ALGORITHM

In this section, an efficient algorithm is proposed to address the mixed optimization problem (6). The purpose of the proposed approach is to reduce the number and the size of QP subproblems to solve in order to reach the optimum. Whereas the brute force strategy previously exposed requires to completely cover the tree of operating mode possibilities represented on figure 1, the basic idea of this algorithm is to prune the tree by cutting the branches that are not able to lead to the optimum.

#### A. Branch cutting

The following notations are used. $u_{dP}^N$ is the discrete valued control subsequence of length $P \le N$ extracted from the discrete valued control sequence $u_{dN}$. The same notation is used for continuous valued ones. The superscript * denotes optimality as for $u_{cN}^*$. However $u_{cP}^{*N}$ is the subsequence of length $P$ extracted from $u_{cN}^*$ and then is non-optimal over a horizon of length $P$.

Consider the following proposition. Given a discrete valued control sequence $u_{dN}$, for any $P \le N$,

$$J_N^*(x(k), u_{dN}) \ge J_P^*(x(k), u_{dP}^N). \quad (10)$$

In the tree of operating mode possibilities (cf. figure 1), on a path from the root to a leaf, cost function (5) still increases.

*Proof:* Split the cost (9) in two terms, one for the cost from time step $k$ to $k + P - 1$ and one for time step $k + P$ to $k + N - 1$:

$$J_N^*(x(k), u_{dN}) = J_P(x(k), (u_{cP}^{*N}, u_{dP}^N)) + \sum_{j=P+1}^{N} L(x(k+j-1), (u_{cP}^{*N}(j), u_{dP}^N(j))). \quad (11)$$

For any $l \in \{k, \dots, k + N - 1\}$, $L(x(l), u(l))$ is positive or null, thus

$$J_N^*(x(k), u_{dN}) \ge J_P(x(k), (u_{cP}^{*N}, u_{dP}^N)). \quad (12)$$

The right hand expression of (12) is the cost given by using the continuous valued control subsequence $u_{cP}^{*N}$ which is extracted from $u_{cN}^*$. This cost is necessarily greater than the optimal cost $J_P^*(x(k), u_{dP}^N)$ obtained by solving a QP subproblem over a horizon of length $P$. Thus the following statement stands:

$$J_P(x(k), (u_{cP}^{*N}, u_{dP}^N)) \ge J_P^*(x(k), u_{dP}^N). \quad (13)$$

Assume to be at a depth $P \le N$ in the tree of possibilities and that the associated cost is greater than a known suboptimum of (5). According to proposition (10) all the remaining branches starting from the considered node can be cut.

#### B. First suboptimum

To improve the efficiency of this branch cutting strategy, it is suitable to begin with a good candidate for the suboptimum. The proposed approach consists in the following basic steps:

* Enumerate at one time step the possibilities for the discrete valued control inputs,

* Solve associated QP subproblems over the prediction horizon according to the current depth in the tree of possibilities,

* Choose the best candidate.

By starting again this procedure, continue the exploration until to reach the end of the prediction horizon, i.e. a suboptimal solution has been found.
<figure><img src="image_2.png" /><figcaption>Fig. 2. Pruning strategy for mixed optimization</figcaption></figure>

### C. Implementation

A tractable approach is to proceed in two steps:

1. Find a suboptimum by using the procedure described above.

2. A suboptimal cost is known. Go one time step backwards to either consider the following best candidate or cut the branches whose cost is already greater than the known suboptimal one.

This pruning strategy is graphically illustrated on figure 2. At bottom, the leaves (symbolized by small circles) are the set of QP sub-problems to solve in order to find the optimum by enumeration (cf. section II-B). The proposed approach consists in going directly towards the suboptimum that is obtained by the series of the optima at one time step enumeration (symbolized by shaded bold circles). Then all the control sequences not yet explored but that are not good candidates (symbolized by empty circles) regarding to the suboptimum are cut (symbolized by crosses) until to find a better suboptimum (symbolized by empty bold circle) or to reach the root of the tree. Only the QP subproblems associated to the nodes marked out by a circle have been needed to be solved. The control sequence obtained that way (symbolized by a square) is the one which minimizes the cost function (5) over the prediction horizon.

## IV. ILLUSTRATIVE EXAMPLE

To illustrate the strategies described in the previous sections, a small size application belonging to the research area of power systems is considered. Voltage control and stability problems are receiving special attention since recently, several major network collapses have occurred around the world, e.g. in North America and Europe in 2003. Power system components contain some interesting mechanisms that makes it suitable for our study. More details about power system components can be found in [9]. Due to the presence of either continuous or discrete valued variables, interesting on voltage stability control issue, power systems

<figure><img src="image_3.png" /><figcaption>Fig. 3. Four-bus power system</figcaption></figure>

belong to the class of SAS.

In the following, the power system example is briefly described and simulation results are given, including efficiency of the used algorithms. All the numerical values are given using the per unit (p.u.) system which is commonly used to deal with this kind of system.

### A. Four-bus power system

Consider the power system shown on figure 3 which is a case study proposed to the control and power engineering community, see [10]. Even though the power system considered here is trivially small from the power engineer's point of view, the main components of a real power system are present:

* $G_1$ is a generator modelling the surrounding network which is assumed to be strong,

* $G_2$ is a generator equipped with a voltage regulator and a field voltage limit,

* $L_1, L_2$ and $L_3$ are transmission lines,

* $C$ is a capacitor bank,

* $T$ is a transformer equipped with a secondary voltage regulator called Load Tap Changer (LTC),

* $l$ is a voltage dependent dynamic load modelling a set of consumers.

Line $L_3$ can be faulted which corresponds to the disconnection of two of its three parallel lines.

The LTC controls the transformer ratio $n_T$ tap by tap within bounds to regulate the secondary voltage. However, since the size of a tap step is quite small (usually in the range of 0.5%-1.5% of the nominal ratio), to simplify the control design, $n_T$ is considered as a continuous variable with $n_T \in [0.9, 1.1]$.

Considered from the viewpoint of power restoration, the load behavior is described by a smooth non-linear
differential equation [11]:

$$ 
\begin{aligned}
T_p \dot{x}_{lp} + x_{lp} &= P_s(V_l) + P_t(V_l) \\
P_l &= (1 - s_l s_{step}) (x_{lp} + P_t(V_l))
\end{aligned} 
\quad (14) $$

where $x_{lp}$ is a continuous state variable, $V_l$ the load voltage magnitude (i.e. $V_l \equiv V_4$), $P_s(V_l) = P_0 V_l^{\alpha_s}$ and $P_t(V_l) = P_0 V_l^{\alpha_t}$ are the steady-state and transient voltage dependencies respectively, $P_l$ is the actual active load power and $T_p$ is the active power recovery time constant. $s_{step}$ is a constant load shedding step and $s_l$ is a discrete control variable. A similar model is used for the reactive load power.

The control inputs are the following:

* $V_4^{ref}$ is the voltage reference of the LTC ($V_4^{ref} \in [0.8, 1.2]$),

* $s_C \in \{0, 1, 2, 3\}$ is the capacitor bank control input (up to 3 parts of 0.1 p.u.),

* $s_l \in \{0, 1\}$ is the load shedding input ($s_l = 1$ corresponds to the disconnection of 5% of the load).

$V_4^{ref}$ is a continuous valued control input whereas $s_C$ and $s_l$ are discrete valued ones. $V_1, V_2, V_3$ and $V_4$ are the voltage magnitude of the corresponding buses and are considered as the outputs.

For this example, the branching factor is equal to 8. To allow some comparisons, a branching factor equal to 4 will be considered too by restricting $s_C \in \{2, 3\}$.

On figure 4, two scenarios are considered to illustrate the system behavior. The first one is the response to a fault that occurs on line $L_3$ at time $t = 150$ s and no control is applied to the system: the voltages collapse. The second one is the response to the same fault when the remaining part of the capacitor bank is connected at time $t = 200$ s as an emergency action: the stability is recovered but the bus voltage $V_3$ level is still unacceptable for power systems operators.

## B. Simulation results

The MPC based method detailed in this paper has been applied to this case study. The control objectives are to keep the voltages close to the nominal value with acceptable levels, i.e. within an interval of $\pm 10\%$ centered around 1 p.u. Power systems are strongly non-linear and in the framework of voltage stability, they can be modeled using a non-linear multimodel approach. A symbolic off-equilibrium linearizing procedure previously developed allows to obtain a SAS model which is valid around the operating point. Time is discretized and using Euler method to approximate the derivative, a discrete time SAS model of power system is so obtained.

<figure><img src="image_4.png" /><figcaption>Fig. 4. Responses to the fault without automatic control</figcaption></figure>

The cost function (5) is used with quadratic terms on $y^T(l) = (V_2(l) V_3(l) V_4(l))^T$ and $u^T(l) = (n_T(l) s_C(l) s_l(l))$. Note that for the sake of simplicity the transformer ratio $n_T$ is used rather than the secondary voltage reference $V_4^{ref}$ in the optimization problem induced by MPC. Thus the implementation of the controller requires an inversion of the transfer between $V_4^{ref}$ and $n_T$ which has been realised by using the prediction of $V_4$ computed by the predictor.

For the control design, the following choices have been made: the time scale of phenomena such as voltage collapse is typically about several minutes and tap changing is subject to the limitation that the transformer ratio can only be changed every 30 s also the sample period has been chosen equal to 30 s; time constants of the load dynamics are about 60 s then the length $N$ of the prediction horizon has been chosen equal to 3. The controller has been tuned by using the following weighing matrices:

$$ Q_y = \begin{pmatrix} 0 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 0.5 \end{pmatrix}, \quad Q_u = \begin{pmatrix} 0 & 0 & 0 \\ 0 & 0 & 0 \\ 0 & 0 & 0.1 \end{pmatrix} \quad (15) $$

Note that due to the presence of bounds on the control inputs, $Q_u$ does not need to be strictly positive defined.

Figure 5 shows that control objectives are achieved: after the occurrence of the fault, the voltage stability is recovered and all the voltages still remain above 0.9 p.u. The remaining part of the capacitor bank has been connected at time $t = 180$ s as shown on figure 6. There has been no need to have resort to load shedding. Furthermore the controller has prevented harmful actions of the LTC on the faulted system by locking it at his current tap which is a well known strategy of power systems operators.
<figure><img src="image_5.png" /><figcaption>Fig. 5. Responses to the fault with MPC</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Fig. 6. Control inputs</figcaption></figure>

### C. Discussion

Brute Force and Mixed Pruning algorithms have been implemented on a common desktop computer running at 2.8 GHz featuring Matlab® R13. The QP solver which is used in the above simulations is quadprog coming from the associated Optimization Toolbox. However powerful commercial solvers are available to deal with quadratic programming. Table I shows average quantities obtained with the four-bus power system example. In the case of the brute force strategy, all the QP subproblems to be solved have the same size that is equal to the dimension of the continuous control input vector multiplied by the length of the prediction horizon: $dim(u_c) \times N$. On the contrary, for the mixed pruning strategy, the size of the QP subproblems starts from $dim(u_c)$ and increases with the depth in the tree of possibilities until to reach $dim(u_c) \times N$ at the first suboptimum (cf. III-B).

<table><caption>TABLE I
ALGORITHM EFFICIENCY</caption><thead><tr><th>Strategy</th><th>Branching Factor</th><th>QP solving</th><th>Computation time (s)</th></tr></thead><tbody><tr><td rowspan="2">Brute Force</td><td>4</td><td>64</td><td>10.8</td></tr><tr><td>8</td><td>512</td><td>77.8</td></tr><tr><td rowspan="2">Pruning</td><td>4</td><td>16</td><td>2.7</td></tr><tr><td>8</td><td>32</td><td>5.3</td></tr></tbody></table>

## V. CONCLUSIONS

The mixed optimization strategies presented in this paper deal with MPC applied to discrete time SAS. As expected, the Brute Force strategy which is simple to implement only suits for (very) small case studies due to the combinatorial explosion and the number of associated QP subproblems that need to be solved in order to reach the optimum. The Mixed Pruning strategy proposed in this paper allows to significantly reduce the size as well as the number of QP subproblems needed to be solved.

## VI. ACKNOWLEDGMENTS

This work has been performed as part of the Energie Unifying Project which deals with technical and economic optimization of power systems (http://www.supelec.fr/ecole/eei/energie).

## REFERENCES

[1] A. Bemporad and M. Morari, “Control of systems integrating logic, dynamics, and constraints,” Automatica, vol. 35, no. 3, pp. 407–427, 1999.

[2] W. Heemels, B. D. Schutter, and A. Bemporad, “On the equivalence of classes of hybrid dynamical models,” in Proceedings of the 40th IEEE Conference on Decision and Control, Orlando, Florida USA, Dec. 2001.

[3] M. Peña, E. Camacho, and S. Piñón, “Hybrid systems for solving model predictive control of piecewise affine system,” in Proceedings of the IFAC Conference on Analysis and Design of Hybrid Systems, Saint-Malo, France, 2003.

[4] R. Alur and G. Pappas, Eds., Hybrid Systems: Computation and Control, ser. Lecture Notes in Computer Science, vol. 2993. Philadelphia, PA, USA: Springer, mar 2004.

[5] B. Kouvaritakis and M. Cannon, Nonlinear predictive control: theory and practice. The Institution of Electrical Engineers, London, 2001.

[6] S. Leirens, J. Buisson, P. Bastard, and J. Coullon, “A hybrid approach to the power system voltage stability issue,” Hybrid System Control Group, Supélec-IETR, Cesson-Sévigné, France, Tech. Rep. ASH04-02, 2004.

[7] G. Goodwin, S. Graebe, and M. Salgado, Control system design. Prentice Hall, 2001.

[8] D. Mayne, J. Rawlings, C. Rao, and P. Scocaert, “Constrained model predictive control: Stability and optimality,” Automatica, vol. 36, no. 6, pp. 789–814, 2000.

[9] P. Kundur, Power System Stability and Control. McGraw-Hill, 1994.

[10] M. Larsson, “The ABB power transmission test case,” Corporate Research ABB Schweiz AG, Tech. Rep., Oct. 2002.

[11] D. Karlsson and D. Hill, “Modelling and identification of nonlinear dynamic loads in power systems,” IEEE Transactions on Power Systems, vol. 9, no. 1, pp. 157–166, Feb. 1994.