# Efficient Model Predictive Algorithms for Tracking of Periodic Signals

Authors: Yun-Chung Chu; Michael Z. Q. Chen
Journal: Journal of Control Science and Engineering
 (2012)

---

## Abstract

This paper studies the design of efficient model predictive controllers for fast-sampling linear time-invariant systems subject to input constraints to track a set of periodic references. The problem is decomposed into a steady-state subproblem that determines the optimal asymptotic operating point and a transient subproblem that drives the given plant to this operating point. While the transient subproblem is a small-sized quadratic program, the steady-state subproblem can easily involve hundreds of variables and constraints. The decomposition allows these two subproblems of very different computational complexities to be solved in parallel with different sampling rates. Moreover, a receding horizon approach is adopted for the steady-state subproblem to spread the optimization over time in an efficient manner, making its solution possible for fast-sampling systems. Besides the conventional formulation based on the control inputs as variables, a parameterization using a dynamic policy on the inputs is introduced, which further reduces the online computational requirements. Both proposed algorithms possess nice convergence properties, which are also verified with computer simulations.


---

## Body

## 1. Introduction

One of the most attractive features of model predictive control (MPC) is its ability to handle constraints [1]. Many other control techniques are conservative in handling constraints, or even try to avoid activating them, thus, sacrificing the best performance that is achievable. MPC, on the contrary, tends to make the closed-loop system operate near its limits and hence produces far better performance. This property of MPC gives it the strength in practice, leading to a wide acceptance by the industry.A very good example of system operating near its limits is a plant being driven by periodic signals to track periodic references. Under this situation, some of the system constraints will be activated repeatedly, and the optimal operating control signal is far from trivial. Just clipping the control signal to fit into the system constraints produces inferior performance typically. And the loss being considered here is not just a transient loss due to sudden disturbances, but indeed a steady-state loss due to a suboptimal operating point. Therefore, the loss is on long term and severe.On the other hand, the successful real-life applications of MPC are mostly on systems with slower dynamics such as industrial and chemical processes [2]. The reason is simply that MPC requires a constrained optimization to be carried out online in a receding horizon fashion [3, 4]. Therefore, to apply MPC to fast-sampling systems, the computational power needed will be substantial. In any case, because of its great success in slow-sampling systems, the trend to extend MPC to fast-sampling systems is inevitable, and many recent researches have been carried out to develop efficient methods to implement MPC in such cases. While some of these works focus on unloading the computational burdens [5–9], others emphasize on code optimization [10–12] and new algorithmic paradigms [13–17].If MPC is applied to the tracking of periodic signals in a receding horizon fashion, the horizon length will be related to the period length, and a long period will imply an online optimization problem of many variables and constraints. For a fast-sampling system, it is essentially an attempt to solve a very big computational problem within a very small time frame. In this paper, we shall analyze the structure of this problem and then propose two efficient algorithms for the task. They aim to make the application of MPC to a fast-sampling system possible by a slight sacrifice on the transient performance, but the optimal or near-optimal steady-state performance of periodic tracking will be maintained.In Section2, the mathematical formation of the problem will be presented. The two algorithms, one based on the concept of receding horizon quadratic programming and the other based on the idea of dynamic MPC policy, will be presented in Sections 3 and 4, respectively. A comment on the actual implementation will be given in Section 5, followed by some computer simulations in Section 6 to illustrate several aspects of the proposed algorithms. Finally, Section 7 concludes the paper.To avoid cumbersome notations likeu(k∣k),u(k+1∣k),…,u(k+Nu-1∣k), the MPC algorithms in this paper will only be presented as if the current time is k=0, and we shall write u(0),u(1),…,u(Nu-1) instead. The reader is asked to bear in mind that the algorithm is actually to be implemented in a receding horizon fashion.


## 2. Problem Formulation

Consider a linear time-invariant plant subject to a periodic disturbance:(1)x+=Ax+B1w+B2u,(2)y=Cx+D1w+D2u,(3)x(0)=x0,
where the superscript + denotes the time-shift operator, that is, (4)x+(k)=x(k+1),
and the disturbance w is measurable and periodic with period Np. The control objective is to construct a control signal u such that the plant output y will track a specific periodic reference r of the same period Np asymptotically with satisfactory transient performance. The control input u is also required to satisfy some linear inequality constraints (e.g., to be within certain bounds). The reference r is not necessarily fixed but may be different for different disturbance w. (For that reason, it may be more appropriate to call w an exogenous signal rather than a disturbance).The algorithms developed in this paper are motivated by the following situations:(1)
the periodNp is very long compared with the typical transient behaviours of the closed-loop system;(2)
the linear inequality constraints onu are persistently active, that is, for any given k̃, there exists a k>k̃ such that u(k) will meet at least one of the associated linear equality constraints;(3)
there is not sufficient computational power to solve the associated quadratic program completely within one sampling interval unless both the control horizon and the prediction horizon are much shorter thanNp.As a matter of fact, without the above considerations and restrictions, the problem is not very challenging and can be tackled with standard MPC approaches for linear systems.The underlying idea of the approach proposed in this paper is that since the transient behaviour of the closed-loop system is expected to be much shorter than the periodNp, we should decompose the original control problem into two: one we call the steady-state subproblem and the other we call the transient subproblem. Hence, the transient subproblem can be solved with a control horizon and a prediction horizon much shorter than Np. While the steady-state subproblem is still very computationally intensive and cannot be solved within one sampling interval, it is not really urgent compared with the transient subproblem, and its computation can be spread over several sampling intervals. Indeed, the two subproblems need not be synchronized even though the transient subproblem depends on the solution to the steady-state subproblem due to the coupled input constraints. The former will utilize the latter whenever the latter is updated and made available to the former. It is only that the transient control will try to make the plant output y track a suboptimal reference if the optimal steady-state control is not available in time.Now let us present the detailed mathematical formulation of our proposed method. Naturally, since bothw and r are periodic with period Np, the solution u should also be periodic of the same period asymptotically, that is, there should exist a periodic signal us(k) such that (5)lim⁡k→∞(u(k)-us(k))=0.
Let xs and ys be the corresponding asymptotic solutions of x and y, and they obviously satisfy the dynamics inherited from (1) and (2):(6)xs+=Axs+B1w+B2us,ys=Cxs+D1w+D2us,xs(0)=xs(Np).
Ideally, we want ys=r but this might not be achievable when us is required to satisfy the specific linear inequality constraints. Therefore, following the spirit of MPC, we shall find us, such that(7)Js=∑kes(k)TQes(k)
is minimized for some positive definite matrix Q, where es(k) is the asymptotic tracking error defined by(8)es=ys-r,
and the summation in (7) is over one period of the signals. This is what we call the steady-state subproblem. In Sections 3 and 4 below, we shall present two different approaches to this steady-state subproblem, with an emphasis on their computational efficiencies.Once the steady-state signals are known, the transient signals defined by(9)ut=u-us,xt=x-xs,yt=y-ys,
satisfy the dynamics(10)xt+=Axt+B2ut,yt=Cxt+D2ut,xt(0)=x0-xs(0),
derived from (1)–(3), subject to the original linear inequality constraints being applied to ut(k)+us(k). Since the control horizon and the prediction horizon for this transient subproblem are allowed to be much shorter than Np, it can be tackled with existing MPC algorithms.It is important to note that in this steady-state/transient decomposition, the steady-state controlus is actually a feedforward control signal determined from w and r, whereas the transient control ut is a feedback control signal depending on x. As an unstable plant can only be stabilized by feedback, but the main interest of the current paper is the computational complexity of the steady-state subproblem, we shall not discuss in depth the stabilization issue, which has been studied quite extensively in the MPC literature. Typically, stabilizability of a constrained system using MPC would be cast as the feasibility of an associated optimization problem. For the numerical example in Section 6, we shall conveniently pick a plant where A is already stable, and hence the following quadratic cost may be adopted for the transient subproblem:(11)Jt=∑k=0Nu-1(yt(k)TQyt(k)+ut(k)TRut(k))+xt(Nu)TPTxt(Nu),
where Nu is the control horizon with Nu≪Np, Q and R are chosen positive definite matrices, and PT is the (weighted) observability gramian obtained from the Lyapunov equation(12)ATPTA-PT+CTQC=0.
The minimization of Jt is simply a standard quadratic program over ut(0),ut(1),…,ut(Nu-1) for a given xt(0). The situation will be more complicated when A is not stable, but one well-known approach is to force the unstable modes to zero at the end of the finite horizon [18].Remark 1.
Essentially, the choice of the cost function (11) with PT from (12) for a stable A means that the projected control action after the finite horizon is set to zero, that is, ut(k)=0 for k≥Nu since the “tail” of the quadratic cost is then given by
(13)∑k=Nu∞yt(k)TQyt(k)=∑k=Nu∞xt(k)TCTQCxt(k)=xt(Nu)TPTxt(Nu).
This terminal policy is valid because the steady-state subproblem has already required that us satisfies the linear inequality constraints imposed on u. Hence, Jt is obviously a Lyapunov function which will be further decreased by the receding horizon implementation when the future control ut(Nu) turns into an optimization variable from zero.Remark 2.
We have deliberately omitted theR-matrix in the steady-state cost Js in (7). The reason is simply that we want to recover the standard linear solution (for perfect asymptotic tracking) as long as us does not hit any constraint.


## 3. Steady-State Subproblem: A Receding Horizon Quadratic Programming Approach

When the periodic disturbancew is perfectly known, the steady-state subproblem is also a conceptually simple (but computationally high-dimensional) quadratic program. One way to know w is simply to monitor and record it over one full period. This, however, does not work well if w is subject to sudden changes. For example, the plant to be considered in our simulations in Section 6 is a power quality device called Unified Power Quality Conditioner [19], where w consists of the supply voltage and the load current of the power system, and both may change abruptly if there are supply voltage sags/swells and variations in the load demands. Indeed, the main motivation of the receding horizon approach in MPC is that things never turn out as expected and the control signal should adapt in an autonomous manner. If the suddenly changed disturbance w can be known precisely only after one full period of observation, the transient response of the steady-state subproblem (not to be confused with the transient subproblem described in Section 2) will be unsatisfactory.One way to overcome this is to introduce an exogenous model for the signalsw and r, as adopted in [19]. Specifically, we construct a state-space model:(14)v+=Avv,(15)w=Cwv,(16)r=Crv,
and assume that both w and r are generated from this model. Since w and r are periodic with period Np, we have(17)AvNp=I.
One simple (but not the only) way to construct Av, as demonstrated similarly in [19] in the continuous-time case, is to make Av a block-diagonal matrix with each block taking the form: (18)[cos⁡(nωTs)-sin⁡(nωTs)sin⁡(nωTs)cos⁡(nωTs)],
where n is an integer, Ts is the sampling time and ωTs×Np=2π. Then the matrices Cw and Cr are just to sum up their respective odd components of v. This essentially performs a Fourier decomposition of the signals w and r, and hence their approximations by Cwv and Crv will be arbitrarily good when more and more harmonics are included in the model.Based on the exogenous model (14)–(17), an observer can be easily constructed to generate (an estimate of) v from the measurements of w and r. From v(0), the model (14)–(17) can then generate predictions of w(0),w(1),…,w(Np-1) and r(0),r(1),…,r(Np-1), and these can be used to find us(0),us(1),…,us(Np-1) by the quadratic program. The use of the exogenous model (14)–(17) typically allows the changed w and r to be identified much sooner than the end of one full period.The quadratic program for the steady-state subproblem can be written as follows:(19)min⁡us(0)Js,Js:=[us(0)v(0)]T[HFFTG][us(0)v(0)],us(0):=[us(0)us(1)⋮us(Np-1)],(20)subjecttoajTus(0)≤bj,j=1,2,…,Nm,
where Nm is the total number of linear inequality constraints. Note that since we assume only input but not state constraints for our original problem, (20) does not depend on v(0) and, hence, the feasibility of any us(0) remains the same even if there is an abrupt change in v(0) (i.e., if v(0) is different from the predicted value from (14) and v(-1)). Furthermore, the active set of constraints remains the same.Definition 3.
The active set𝒜(us(0)) of any feasible us(0) satisfying (20) is the subset of {1,2,…,Nm} such that j∈𝒜(us(0)) if and only if ajTus(0)=bj.Next, we present a one-step active set algorithm to solve the quadratic program (19)-(20) partially.Algorithm 4.
Given an initial feasibleus(0) and a working set 𝒲0⊂𝒜(us(0)).Let the set of working constraints
(21)ajTus(0)≤bj,j∈W0,
be represented by
(22)A0us(0)≤b0,
where the inequality sign applies componentwise, that is, each row of A0, b0 represents a working constraint in (21). (1)
Compute the gradient(23)g0=Hus(0)+Fv(0),
and the null space of A0, denoted Z0 by
(24)A0Z0=0.
If Z0Tg0≈0, go to step  (5).(2)
Compute a search directionw0=Z0ŵ0, where
(25)(Z0THZ0)ŵ0+Z0Tg0=0.(3)
Let(26)α0=min⁡j∈A(us(0))∖W0s.t.ajTw0>0.⁡bj-ajTus(0)ajTw0.(4)
Ifα0≥1, go to step  (5). Otherwise, update us(0) to us*(0) by
(27)us*(0)=us(0)+α0w0,
and add a (blocking) constraint to 𝒲0 to form a new working set 𝒲0* according to the method described in Remark 7 below. Quit.(5)
Updateus(0) to us*(0) by
(28)us*(0)=us(0)+w0.
Compute the Lagrange multiplier λ0 from
(29)A0Tλ0+g0=0
to see whether any component of λ0 is negative. If yes, remove one of the constraints corresponding to a negative component of λ0 from 𝒲0 to form a new working set 𝒲0* according to the method described in Remark 7 below. Quit.Algorithm4 can be interpreted as follows. It solves the equality-constrained quadratic program:(30)min⁡us′(0)Js′,Js′:=[us′(0)v(0)]T[HFFTG][us′(0)v(0)]subjecttoajTus′(0)=bj,j∈W0,,
and then searches along the direction(31)w0=us′(0)-us(0),
until it is either blocked by a constraint not in 𝒲0 (step  (4)) or the optimal us′(0) is reached (step  (5)). This is indeed a single step of the standard active set method (the null-space approach) for quadratic programming [20, 21] except for the modifications that will be detailed in Remark 7 below. In other words, if we apply Algorithm 4 repeatedly to the new us*(0), 𝒲*(0), it will converge to the solution of the original inequality-constrained quadratic program (19)-(20) within a finite number of steps, and the cost function Js is strictly decreasing. However, here we only apply a single step of the active set method due to the limited computational power available within one sampling interval. Furthermore, we do not even assume that the single step of computation can be completed within Ts. Let NaTs be the time required or allowed to carry out Algorithm 4. To complete the original quadratic program (19)-(20) in a receding horizon fashion, we need to forward us*(0), 𝒲*(0) to us(Na), 𝒲(Na) by rotating the components of us*(0) by an amount of Na since it is supposed to be periodic:(32)[us(Na)us(Na+1)⋮us(Np-1)us(Np)us(Np+1)⋮us(Na+Np-1)]=[us*(Na)us*(Na+1)⋮us*(Np-1)us*(0)us*(1)⋮us*(Na-1)].
Obviously, Algorithm 4 will then continue to solve an equivalent quadratic program as long as v strictly follows the exogenous dynamics (14). Hence we have the following convergence result.Proposition 5.
Algorithm4 together with the rotation of the components of us*(0) in (32) will solve the quadratic program (19)-(20) in finite time as long as v(k) satisfies the exogenous dynamics (14).Proof.
From the argument above it is easy to see that as long asv(k) follows the dynamics (14), the algorithm is consistently solving essentially the same quadratic program. So it remains to check that the convergence proof of the standard active set algorithm remains valid despite the modifications we shall detail in Remark 7, which is indeed the case.Of course, the most interesting feature of the receding horizon approach is that the solution will adapt autonomously to the new quadratic program if there is an abrupt change inv. Since constraint (20) is independent of v, an abrupt change in v will not destroy the feasibility of the previously determined us and the working set 𝒲 determined previously also remains a valid subset of the active set. Hence, the receding horizon active set method will continue to work even though the objective function (19) has changed. However, if it is necessary to include not only control but also state constraints into the original problem formulation, we shall then require a quadratic programming algorithm (other than the active set method in its simplest form) that does not assume the initial guess to be feasible.Remark 6.
There could be two possible ways to arrange the steps in Algorithm4. One is to update the working set 𝒲 followed by us, and the other is to update us followed by the working set 𝒲. In the receding horizon framework, it might seem at first glance that the first choice is the right one, since we shall then avoid basing the optimization of us on an “outdated” working set if v happens to have changed. However, it turns out that the first choice is actually undesirable. One well-known “weakness” of the active set method is that it is not so easy to remove a constraint once it enters the working set 𝒲. This becomes even more a concern in the receding horizon framework. If v has changed, and so has the objective function (19), the original stationary us* obtained in step  (5) may no longer be stationary, and it will require at least an additional iteration to identify the new stationary point before we can decide whether any constraint can be dropped from the working set or not. This will seriously slow down the transient response of the steady-state subproblem. Indeed, once v has changed, many of the constraints in the previous working set are no longer sensible, and it will be wiser to drop them hastily rather than being too cautious only to find much later that the constraints should still be dropped eventually.Remark 7.
One key element in the active set method of the quadratic program is to add or drop a constraint to or from the working set𝒲. The constraint being added belongs to the blocking set ℬ, defined as those constraints corresponding to the minimum α0 in (26). Physically, they are the constraints that will be first encountered when we try to move us(0) to us′(0) in (31). The constraint being dropped belongs to the set ℒ, defined as those constraints corresponding to a negative component of the Lagrange multiplier in (29). The active set method will converge in finite time no matter which constraint in ℬ will be added or which constraint in ℒ will be dropped. One standard and popular choice in the conventional active set method is that the one in ℒ corresponding to the most negative component of λ will be dropped, whereas the choice from ℬ will be arbitrary. This is a very natural choice when there is no other consideration.
However, in the receding horizon framework, one other (and in fact important) consideration emerges, which is the execution time of the control input(s) associated with a constraint. Specifically, if Algorithm4 takes time NaTs to carry out, then 𝒲0* updated in the current iteration will be used to optimize us(Na) in the next iteration,
(33)us(Na):=[us(Na)us(Na+1)⋮us(Na+Np-1)],
but the outcome of that optimization is ready only at k=2Na, based on which the transient control ut is computed. Suppose that the transient subproblem takes one sampling interval to solve, then the transient subproblem at k=2Na will update u(2Na+1)=us(2Na+1)+ut(2Na+1) (see Section 5 below for a more detailed discussion). Hence, the “time priority” of us is 2Na+1,2Na+2,…,Np-1,0,1,…,2Na and from this argument, we choose to drop the constraint in ℒ that is associated with the first us in this sequence or to add the constraint in ℬ that is associated with the last us in this sequence (of course the most negative Lagrange multiplier can still be used as the second criterion if two constraints in ℒ happen to have the same urgency). The proposal here aims to assign most freedom to the most urgent control input in the optimization, which makes sense in the receding horizon framework since the less urgent inputs may be reoptimized later.Remark 8.
Basically, the approach proposed in this section is to spread the original quadratic program over many intervals, so that each interval only carries out one iteration of the algorithm, and also to ensure that the quadratic program being solved is consistent when the prediction of the exogenous model is valid, but will migrate to a new quadratic program when there is a change inv(k). It is worth mentioning that the original standard MPC is a static controller by nature, since the true solution of a complete quadratic program is independent of the MPC’s decisions in the past (past decisions can help to speed up the computations but will not affect the solution), but by spreading the algorithm over several intervals, it is turned into a dynamic controller with internal state us(k),𝒲(k).


## 4. Steady-State Subproblem: A Dynamic Policy Approach

The approach proposed in Section3 optimizes us. Consequently, the number of (scalar) variables being optimized is proportional to Np. To further cut down the computations required, this section proposes another approach based on the idea of a dynamic policy, inspired by the works of [13, 22, 23]. This approach optimizes a smaller number of variables typically, and the number is independent from Np, although the number of linear inequality constraints remains the same. In return, the optimization result is expected to be slightly inferior to that of Section 3 due to the reduction of variables (degree of freedom). However, it should be noted that the number of optimized variables in this second approach is adjustable based upon the designer’s wish.The central idea of the dynamic policy approach [13, 22, 23] is that instead of optimizing the control input directly, we generate the control input by a dynamic system of which the initial system state is optimized. This is similar to what we have done to w and r in Section 3. Specifically, we assume us is also generated from a state-space model:
(34)v̂+=Av̂v̂,(35)us=Cv̂v̂,
where(36)Av̂Np=I.
This state-space model is designed a priori but the initial state v̂(0) will be optimized online. Obviously, the quadratic program (19)-(20) becomes(37)min⁡v̂(0)Ĵs,Ĵs:=[v̂(0)v(0)]T[ĤF̂F̂TG][v̂(0)v(0)],(38)subjecttoâjTv̂(0)≤bj,j=1,2,…,Nm,
where(39)Ĥ=ÔTHÔ,F̂=FÔ,âj=ajÔ,j=1,2,…,Nm,
and 𝒪̂ is the observability matrix(40)Ô:=[Cv̂Cv̂Av̂⋮Cv̂Av̂Np-1].
The number of variables in this new quadratic program is the dimension of v̂(0), denoted by nv̂. If Av̂ is constructed from the method of Fourier decomposition described in Section 3, Shannon’s sampling theorem implies that a sufficiently large but finite nv̂ will guarantee a full reconstruction of the original optimization variable us(0). On the other hand, a smaller nv̂ restricts the search to a lower dimensional subspace of us(0) and hence the optimization is easier but suboptimal.One natural choice of the dynamicsAv̂ is to make Av̂=Av in the exogenous model (14)–(17). Of course, it should be noted that constrained control is generally a nonlinear problem, and therefore the number of harmonics to be included in us may exceed that of w and r in order to achieve the true optimal performance. However, we could have over-designed the exogenous model (14)–(17) to include more harmonics in Av than necessary for w and r, thus making the choice Av̂=Av here not so conservative. The simulation results in Section 6 will demonstrate both cases.It remains to choose the matrixCv̂ in (35). The one we suggest here is based on the linear servomechanism theory [24–26], which solves the linear version of our problem when there is no input constraint. Essentially, when there is no constraint, perfect asymptotic tracking (i.e., ys=r or es=0) can be achieved by solving the regulator equation:(41)XAv=AX+B1Cw+B2U,Cr=CX+D1Cw+D2U,
for the matrices X, U, and then let(42)us=Uv,
which also implies(43)xs=Xv.
Therefore, to recover the optimal (or perfect) solution in the linear case when us does not hit any constraint, the state-space model of us may be chosen as(44)ṽ+=Avṽ,us=Uṽ.
However, this state-space model is not guaranteed to be observable. When it is not, the resulting Ĥ in (37) will become semidefinite instead of strictly positive definite. To overcome this, we suggest to perform an orthogonal state transformation(45)[*v̂]=Tṽ,
to bring (44) to the Kalman decomposition form(46)[*v̂]+=[**0Av̂][*v̂],us=[0Cv̂][*v̂],
and hence obtain a reduced-order model to become (34)-(35). It is easy to verify that Av̂Np=I since (47)[**0Av̂]
is upper block-triangular and (48)[**0Av̂]Np=(TAvT-1)Np=TAvNpT-1=I.Certainly, the discussion above only provides a suggestion of how to choose the state-space model forus, which we shall also adopt for our simulations in Section 6, but, in general, the designer should feel free to employ any valid state-space model to suit his problem.Remark 9.
Having reparameterized the quadratic program in terms ofv̂(0) rather than us(0), we can apply a similar version of Algorithm 4  to (37)-(38). In other words, it is not necessary to solve the quadratic program completely within one sampling interval. Instead of rotating the components of us*(0) to obtain us(Na), we obtain v̂(Na) by Av̂Nav̂*(0).


## 5. Implementation Issues and Impacts on Transient Performance

Before we present the simulation results, let us comment on the impact of computational delay on the transient subproblem in Section2. First of all, we assume that the transient quadratic program can be solved completely within one sampling interval. Therefore, despite the way we presented the cost function Jt in (11), in actual implementation we shall optimize ut(1),ut(2),…,ut(Nu) based on xt(1) at time k=0, instead of ut(0),ut(1),…,ut(Nu-1) based on xt(0). The unknown xt(1) can be projected from the known variables and system dynamics. After the optimization is carried out to obtain ut(1),ut(2),…,ut(Nu), the control input to be executed at k=1 will be u(1)=us(1)+ut(1). Bear in mind that all these calculations can only be based on the best knowledge of the signals at k=0.Figure1 summarizes how the unknown variables can be computed from the known variables. The variables in each layer are derived from those variables directly below them, but in actual implementation it is sometimes possible to derive explicit formulas to compute the upper layer from the lower layer bypassing the intermediate layer, thus not requiring those intermediate calculations online. The variable on top is the control action u(1), computed from the variables of the steady-state subproblem on the left, and those of the transient subproblem on the right, separated by a solid line. The variables in the bottom layer, v(0) is provided by the observer described in Section 3, x(0) is a measurement of the current plant state, and u(0) is the control input calculated by the algorithm at k=-1. Note that to compute ut(1), the values of us(1),us(2),…,us(Nu) are needed to form the constraints for the transient quadratic program since the original linear inequality constraints apply to u(k)=us(k)+ut(k). On the other hand, xs(1) can be written as a linear function of v(1) and us(1) (or v̂(1)) explicitly. Finally, the steady-state subproblem requires a computational time of NaTs, implying that the solution us(0) provided by the steady-state subproblem at k=0 is based on a measurement of v(k) at some k between -2Na+1 and -Na. So in the worst case, u(1) is based on some information as old as v(-2Na+1), which corresponds to a worst-case delay of 2NaTs. For instance, if Na=1, the control u(k) is computed from a measurement of x(k-1), v(k-1), and v(k-2).Figure 1
Derivations of unknown variables from known variables.Remark 10.
Although we said in Section2 that the transient subproblem was not the main interest of this paper, it is an ideal vehicle to demonstrate the power of MPC since the “useful freedom” of ut(k) may have been totally consumed by us(t) when the latter hits a constraint. For example, the simulations to be discussed in Section 6 have the constraint
(49)|u(k)|=|us(k)+ut(k)|≤1.
If us(k) already saturates at ±1, one side of ut(k) is lost but that could be the only side to make the reduction of Jt possible, thus forcing ut(k) to zero. So the input constraint does not only restrict the magnitude, but also the time instant to apply the transient control ut. Such problem is extremely difficult for conventional control techniques.


## 6. Simulation Results

In this section we use an example to demonstrate the performance of our algorithms by computer simulations. The plant is borrowed from [19] and represents a power quality device called Unified Power Quality Conditioner (UPQC), which has the following continuous-time state-space model: (50)ẋ=Ax+B1w+B2u,y=Cx+D1w+D2u,A=[-RlLl00-1Ll-1Ll0-RseLse0-1Lse000-RshLsh0-1Lsh1Cse1Cse0001Csh01Csh00],B1=[1Ll00000000-1Csh],B2=[00Vdc2Lse00Vdc2Lsh0000],C=[0000110000],D1=[0000],D2=[0000].
The exogenous input w is composed of the supply voltage and the load current, which are periodic at 50 Hz but may consist of higher-order harmonics. The plant output y is composed of the load voltage and the supply current, which will be made to track designated pure sine waves of 50 Hz. The control input u is composed of two switching signals across the voltage source inverters (VSIs), both of which are required to satisfy the bounds -1≤u≤1. The general control objective is to maintain y to the desired waveforms despite possible fluctuations in w like supply voltage sags/swells or load demand changes. To apply the MPC algorithms proposed, we obtain a discrete-time version of the above state-space model by applying a sampling interval of Ts=0.2ms (i.e., 100 samples per period). Small-sized quadratic programs (such as our transient subproblem) can possibly be solved within such a short time thanks to the state-of-the-art code optimization [12], which reports sampling rates in the range of kHz, but to solve a big quadratic program like our steady-state subproblem we shall resort to the technique of Algorithm 4. Note that in our formulation, the transient subproblem and the steady-state subproblem can be solved in parallel. Although the optimization of ut depends on us, the transient control ut(k+1) is computed from us(k) which is made available by the steady-state subproblem in the previous step. So it is independent of the current steady-state subproblem being solved.As typical in a power system, we assume only odd harmonics in the supply voltage and the load current. Hence we can reduce the computations in the steady-state subproblem by the following easy modifications from the standard algorithms presented in Sections3 and 4; Np may be chosen to represent half of the period instead of the whole period, satisfying(51)xs(0)=-xs(Np),
instead of (6), and(52)AvNp=-I,Av̂Np=-I,
instead of (17) and (36), with ωTs×Np=π instead of 2π. The rotation operation in (32) should also become(53)[us(Na)us(Na+1)⋮us(Np-1)us(Np)us(Np+1)⋮us(Na+Np-1)]=[us*(Na)us*(Na+1)⋮us*(Np-1)-us*(0)-us*(1)⋮-us*(Na-1)].
This cuts down half of the scalar variables as well as constraints in the quadratic program.The model parameters of the UPQC used in our simulations are summarized in Table1 for the circuit components and Table 2 for the line and VSI impedances. They are the same as those values in [19], except for Vdc which we have changed from 400 V to 320 V so as to produce a saturated control u more easily. Note that Vdc is the DC-link voltage, which determines how big a fluctuation in the supply voltage or load current the UPQC can handle. In other words, saturation occurs when the UPQC is trying to deal with an unexpected voltage sag/swell or load demand that is beyond its designed capability.Table 1
Values of the components of the UPQC.
ComponentVdcLseLshCseCshValue320 V5.0 mH1.2 mH10μF20μFTable 2
Line impedance and VSI impedances of the UPQC.
ComponentRlLlRseRshValue0.01Ω1.0 mH0.01Ω0.01ΩThe simulation scenario is summarized in Figure2. Both the supply voltage and the load current consist of odd harmonics up to the 9th order. Despite the harmonics, it is desirable to regulate the load voltage to a fixed pure sine wave, whereas the supply current should also be purely sinusoidal, but its magnitude and phase are selected to maintain a power factor of unity and to match the supply active power to the active power demanded by the load, which means the reference of this supply current is w-dependent. The waveforms of both w and r are shown in Figure 2. The simulation scenario is designed such that the steady-state control us is not saturated at the beginning. At t=0.5s, a voltage sag occurs which reduces the supply voltage to 10% of its original value. The UPQC is expected to keep the load voltage unchanged but (gradually) increase the supply current so as to retain the original active power. This will drive us into a slightlysaturated situation. At t=1.0s, the load demand increases, causing the reference of the supply current to increase again, and us will become deeplysaturated. At t=1.5s, the voltage sag is cleared and the supply voltage returns to its starting value, but the (new) load demand remains. Although the load demand is still higher than its initial value, the us required will just be within the bounds of ±1, thus leaving the saturation region to return to the linear region. So, in short, us is expected to experience “linear → slightly saturated →  deeply saturated →linear but nearly saturated” in this simulation scenario.Figure 2
Simulation scenario. Voltage sag att=0.5s; load demand changed at t=1.0s; sag cleared at t=1.5s.To evaluate the performance of our algorithms, we compare them to two other cases. In the first case, instead of Algorithm4, the complete quadratic program is solved in each iteration of the steady-state subproblem every Na-sampling intervals. We call this case the complete QP, and it serves to indicate how much transient performance has been sacrificed (in theory) by spreading the quadratic program over a number of iterations. In the second case, the constraints are totally ignored, such that the optimal us(0) in the steady-state subproblem is supposed to be(54)us(0)=[UUAv⋮UAvNp-1]v(0),
where U is the solution to the regulator equation (41). The transient subproblem can also be solved by(55)ut=Fxt,
where F is the optimal state-feedback gain minimizing the transient quadratic cost (56)∑k=0∞(yt(k)TQyt(k)+ut(k)TRut(k)).
However, the combined input u is still clipped at ±1. We label this control law the multivariable regulator (MVR) following the linear servomechanism theory. This case serves to indicate how bad the quadratic cost Js can be if a linear control law is used without taking the constraints into consideration.Note that in both cases, the computational delays discussed in Section5 will be in force, where u(k+1) instead of u(k) will be optimized and the steady-state control us is only updated every Na-sampling intervals. In reality, of course the MVR should involve negligible computational delay, whereas the complete QP should need a longer time to solve than Algorithm 4, but we are merely using their associated quadratic costs here to analyze the behaviours of our algorithms.Figure3 plots the steady-state cost Js of our first approach in Section 3 based on receding horizon quadratic programming together with the costs in the other two cases. Na is assumed to be 3 in this simulation. The transient subproblem has a control horizon of Nu=5, corresponding to 10 scalar variables and 20 scalar inequality constraints. On the other hand, the steady-state subproblem has Np=50 corresponding to 100 variables and 200 constraints. As shown in Figure 3, all Js are zero prior to t=0.5s and should also settle down to zero after t=1.5s. It is observed that the transient response of our Js is pretty close to that of the complete QP during the transitions from “linear” to “slightly saturated" and from “slightly saturated" to “deeply saturated,” but is poorer when it tries to return from “deeply saturated” to “linear.” This can probably be attributed to the weakness of the active set method in removing a constraint from the working set as discussed in Remark 6. Figure 3 also indicates that the MVR settles down to a much higher Js value when a saturation occurs, due to its ignorance of the control constraints. The exact values of Js just prior to t=1.0s and t=1.5s are summarized in Table 3, which are about 6-7 times the Js values of our algorithm.Table 3
Summary ofJs values for various cases studied.
Steady-state costJsNaNumber ofatt<1.0satt<1.5sVariablesConstraintsMVR0.115653.276MPC, receding horizon quadratic programming0.01788.75273100 (Np×2)200 (Np×2×2)MPC, dynamic policy (up to 9th harmonics)0.02259.02731nv̂=20 (5×2×2)200 (Np×2×2)MPC, dynamic policy (up to 29th harmonics)0.01828.76402nv̂=60 (15×2×2)200 (Np×2×2)Transient subproblem10 (Nu×2)20 (Nu×2×2)Figure 3
Quadratic costJs of the steady-state subproblem, the receding horizon quadratic programming approach.Figure4 zooms into the first component of the control input u. The top two plots draw attentions to the steady-state control us (u is essentially us just prior to t=1.0s and t=1.5s). It is observed that the differences in us between MVR and MPC are very subtle. Compared with the MVR, the MPC us merely reaches and leaves its limit (+1 here) at very slightly different time instants and also produces some “optimized ripples” of less than 0.25% around that limit instead of a “flat” value as adopted in the clipped linear control law, but by doing these little things the MPC manages to bring Js down by almost one order of magnitude. This demonstrates how nontrivial the optimal us can be. We can also see from the plots that only one constraint is active in the “slightly saturated” situation whereas multiple constraints are active in the “deeply saturated” saturation. On the other hand, the bottom two plots in Figure 4 illustrates our discussion in Remark 10. The plots clearly show that during certain moments of the transient stages (t>1.0s and t>1.5s), the transient control ut is “disabled” due to the saturation of the steady-state control us. Note that we are labeling the dashed blue curve as “us used to compute ut” since it is slightly different from the actual us. For instance, ut(k+1) is computed from the knowledge of us at time k, which is not exactly the same as us(k+1). Obviously, ut is not just disabled whenever us saturates. It happens only when the desired direction of ut violates the active constraint.Figure 4
The first component of controlsu and ut before and after transitions at t=1.0s and t=1.5s.Next, let us look at the performance of our second MPC approach in Section4 based on dynamic policy. Odd harmonics up to the 9th order are included in Av̂ resulting in a total of nv̂=20 variables. See Table 3. Since the number of variables is much lower than that of the first approach, we assume Na=1 here, that is, one iteration of Algorithm 4 (equivalent version) is carried out in each sampling interval, whereas the transient subproblem is solved completely within each sampling interval. The transient performance of Js is plotted in Figure 5. Note that the MVR curve exhibits a slightly different transient from Figure 3 since their Na values are different. The dynamic policy approach clearly shows a faster transient response than the receding horizon quadratic programming approach, not only because of a smaller Na but also a smaller-sized quadratic program overall. However, the drawback is a slightly suboptimal Js, as indicated in Table 3.Figure 5
Quadratic costJsof the steady-state subproblem, the dynamic policy approach.As mentioned in Section4, it is possible to over-design Av, and hence Av̂, so that the optimal Js in this second MPC method will approach the first MPC method. For example, although we only have odd harmonics up to the 9th order in w, we may include odd harmonics up to the 29th order in Av and Av̂. The results are also recorded in Table 3, and we see that this Js value is very close to the optimal one in the first method.


## 7. Conclusions

To apply MPC to fast-sampling systems with input constraints for the tracking of periodic references, efficient algorithms to reduce online computational burdens are necessary. We have decomposed the tracking problem into a computationally complex steady-state subproblem and a computationally simple transient subproblem, and then proposed two approaches to solve the former. The first approach, based on the concept of receding horizon quadratic programming, spreads the optimization over several sampling intervals, thus reducing the computational burdens at the price of a slower transient response. The second approach, based on the idea of a dynamic policy on the control input, further reduces the online computations at the price of a slightly suboptimal asymptotic performance. Despite the limitations, these approaches make the application of MPC to fast-sampling systems possible. Their transient behaviours and steady-state optimality have been analyzed via computer simulations, which have also demonstrated that the steady-state subproblem and the transient subproblem can be solved in parallel with different sampling rates. When the methods proposed in this paper are combined with modern code optimizations, the applicability of MPC to the servomechanism of fast-sampling constrained systems will be greatly enhanced.



---

*Source: 729748-2012-01-18.xml*