Journal of Control Science and Engineering

Journal of Control Science and Engineering / 2012 / Article
Special Issue

Model Predictive Control

View this Special Issue

Research Article | Open Access

Volume 2012 |Article ID 729748 |

Yun-Chung Chu, Michael Z. Q. Chen, "Efficient Model Predictive Algorithms for Tracking of Periodic Signals", Journal of Control Science and Engineering, vol. 2012, Article ID 729748, 13 pages, 2012.

Efficient Model Predictive Algorithms for Tracking of Periodic Signals

Academic Editor: Baocang Ding
Received30 Jun 2011
Accepted31 Aug 2011
Published18 Jan 2012


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.

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 Section 2, 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 like 𝑢(𝑘∣𝑘),𝑢(𝑘+1∣𝑘),…,𝑢(𝑘+𝑁𝑢−1∣𝑘), the MPC algorithms in this paper will only be presented as if the current time is 𝑘=0, and we shall write 𝑢(0),𝑢(1),…,𝑢(𝑁𝑢−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𝑤+𝐵2𝑢,(1)𝑦=𝐶𝑥+𝐷1𝑤+𝐷2𝑢,(2)𝑥(0)=𝑥0,(3) where the superscript + denotes the time-shift operator, that is, 𝑥+(𝑘)=𝑥(𝑘+1),(4) and the disturbance 𝑤 is measurable and periodic with period 𝑁𝑝. The control objective is to construct a control signal 𝑢 such that the plant output 𝑦 will track a specific periodic reference 𝑟 of the same period 𝑁𝑝 asymptotically with satisfactory transient performance. The control input 𝑢 is also required to satisfy some linear inequality constraints (e.g., to be within certain bounds). The reference 𝑟 is not necessarily fixed but may be different for different disturbance 𝑤. (For that reason, it may be more appropriate to call 𝑤 an exogenous signal rather than a disturbance).

The algorithms developed in this paper are motivated by the following situations: (1)the period 𝑁𝑝 is very long compared with the typical transient behaviours of the closed-loop system; (2)the linear inequality constraints on 𝑢 are persistently active, that is, for any given ̃𝑘, there exists a ̃𝑘𝑘> such that 𝑢(𝑘) 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 than 𝑁𝑝.

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 period 𝑁𝑝, 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 𝑁𝑝. 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 𝑦 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 both 𝑤 and 𝑟 are periodic with period 𝑁𝑝, the solution 𝑢 should also be periodic of the same period asymptotically, that is, there should exist a periodic signal 𝑢𝑠(𝑘) such that limğ‘˜â†’âˆžî€·ğ‘¢(𝑘)−𝑢𝑠(𝑘)=0.(5) Let 𝑥𝑠 and 𝑦𝑠 be the corresponding asymptotic solutions of 𝑥 and 𝑦, and they obviously satisfy the dynamics inherited from (1) and (2):𝑥+𝑠=𝐴𝑥𝑠+𝐵1𝑤+𝐵2𝑢𝑠,𝑦𝑠=𝐶𝑥𝑠+𝐷1𝑤+𝐷2𝑢𝑠,𝑥𝑠(0)=𝑥𝑠𝑁𝑝.(6) Ideally, we want 𝑦𝑠=𝑟 but this might not be achievable when 𝑢𝑠 is required to satisfy the specific linear inequality constraints. Therefore, following the spirit of MPC, we shall find 𝑢𝑠, such that𝐽𝑠=𝑘𝑒𝑠(𝑘)T𝑄𝑒𝑠(𝑘)(7) is minimized for some positive definite matrix 𝑄, where 𝑒𝑠(𝑘) is the asymptotic tracking error defined by𝑒𝑠=𝑦𝑠−𝑟,(8) 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) satisfy the dynamics𝑥+𝑡=𝐴𝑥𝑡+𝐵2𝑢𝑡,𝑦𝑡=𝐶𝑥𝑡+𝐷2𝑢𝑡,𝑥𝑡(0)=𝑥0−𝑥𝑠(0),(10) derived from (1)–(3), subject to the original linear inequality constraints being applied to 𝑢𝑡(𝑘)+𝑢𝑠(𝑘). Since the control horizon and the prediction horizon for this transient subproblem are allowed to be much shorter than 𝑁𝑝, it can be tackled with existing MPC algorithms.

It is important to note that in this steady-state/transient decomposition, the steady-state control 𝑢𝑠 is actually a feedforward control signal determined from 𝑤 and 𝑟, whereas the transient control 𝑢𝑡 is a feedback control signal depending on 𝑥. 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 𝐴 is already stable, and hence the following quadratic cost may be adopted for the transient subproblem:𝐽𝑡=𝑁𝑢−1𝑘=0𝑦𝑡(𝑘)T𝑄𝑦𝑡(𝑘)+𝑢𝑡(𝑘)T𝑅𝑢𝑡(𝑘)+𝑥𝑡𝑁𝑢T𝑃𝑇𝑥𝑡𝑁𝑢,(11) where 𝑁𝑢 is the control horizon with 𝑁𝑢≪𝑁𝑝, 𝑄 and 𝑅 are chosen positive definite matrices, and 𝑃𝑇 is the (weighted) observability gramian obtained from the Lyapunov equation𝐴T𝑃𝑇𝐴−𝑃𝑇+𝐶T𝑄𝐶=0.(12) The minimization of 𝐽𝑡 is simply a standard quadratic program over 𝑢𝑡(0),𝑢𝑡(1),…,𝑢𝑡(𝑁𝑢−1) for a given 𝑥𝑡(0). The situation will be more complicated when 𝐴 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 𝑃𝑇 from (12) for a stable 𝐴 means that the projected control action after the finite horizon is set to zero, that is, 𝑢𝑡(𝑘)=0 for 𝑘≥𝑁𝑢 since the “tail” of the quadratic cost is then given by âˆžî“ğ‘˜=𝑁𝑢𝑦𝑡(𝑘)T𝑄𝑦𝑡(𝑘)=âˆžî“ğ‘˜=𝑁𝑢𝑥𝑡(𝑘)T𝐶T𝑄𝐶𝑥𝑡(𝑘)=𝑥𝑡𝑁𝑢T𝑃𝑇𝑥𝑡𝑁𝑢.(13) This terminal policy is valid because the steady-state subproblem has already required that 𝑢𝑠 satisfies the linear inequality constraints imposed on 𝑢. Hence, 𝐽𝑡 is obviously a Lyapunov function which will be further decreased by the receding horizon implementation when the future control 𝑢𝑡(𝑁𝑢) turns into an optimization variable from zero.

Remark 2. We have deliberately omitted the 𝑅-matrix in the steady-state cost 𝐽𝑠 in (7). The reason is simply that we want to recover the standard linear solution (for perfect asymptotic tracking) as long as 𝑢𝑠 does not hit any constraint.

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

When the periodic disturbance 𝑤 is perfectly known, the steady-state subproblem is also a conceptually simple (but computationally high-dimensional) quadratic program. One way to know 𝑤 is simply to monitor and record it over one full period. This, however, does not work well if 𝑤 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 𝑤 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 𝑤 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 signals 𝑤 and 𝑟, as adopted in [19]. Specifically, we construct a state-space model:𝑣+=𝐴𝑣𝑣,(14)𝑤=𝐶𝑤𝑣,(15)𝑟=𝐶𝑟𝑣,(16) and assume that both 𝑤 and 𝑟 are generated from this model. Since 𝑤 and 𝑟 are periodic with period 𝑁𝑝, we have𝐴𝑁𝑝𝑣=𝐼.(17) One simple (but not the only) way to construct 𝐴𝑣, as demonstrated similarly in [19] in the continuous-time case, is to make 𝐴𝑣 a block-diagonal matrix with each block taking the form: cos𝑛𝜔𝑇𝑠−sin𝑛𝜔𝑇𝑠sin𝑛𝜔𝑇𝑠cos𝑛𝜔𝑇𝑠,(18) where 𝑛 is an integer, 𝑇𝑠 is the sampling time and 𝜔𝑇𝑠×𝑁𝑝=2𝜋. Then the matrices 𝐶𝑤 and 𝐶𝑟 are just to sum up their respective odd components of 𝑣. This essentially performs a Fourier decomposition of the signals 𝑤 and 𝑟, and hence their approximations by 𝐶𝑤𝑣 and 𝐶𝑟𝑣 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) 𝑣 from the measurements of 𝑤 and 𝑟. From 𝑣(0), the model (14)–(17) can then generate predictions of 𝑤(0),𝑤(1),…,𝑤(𝑁𝑝−1) and 𝑟(0),𝑟(1),…,𝑟(𝑁𝑝−1), and these can be used to find 𝑢𝑠(0),𝑢𝑠(1),…,𝑢𝑠(𝑁𝑝−1) by the quadratic program. The use of the exogenous model (14)–(17) typically allows the changed 𝑤 and 𝑟 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:min𝐮𝑠(0)𝐽𝑠,𝐽𝑠𝐮∶=𝑠(0)𝑣(0)T𝐅𝐇𝐅T𝐆𝐮𝑠,𝐮(0)𝑣(0)ğ‘ âŽ¡âŽ¢âŽ¢âŽ£ğ‘¢(0)∶=𝑠𝑢(0)𝑠⋮𝑢(1)ğ‘ î€·ğ‘ğ‘î€¸âŽ¤âŽ¥âŽ¥âŽ¦,−1(19)subjectto𝐚T𝑗𝐮𝑠(0)≤𝐛𝑗,𝑗=1,2,…,𝑁𝑚,(20) where 𝑁𝑚 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 𝑣(0) and, hence, the feasibility of any 𝐮𝑠(0) remains the same even if there is an abrupt change in 𝑣(0) (i.e., if 𝑣(0) is different from the predicted value from (14) and 𝑣(−1)). Furthermore, the active set of constraints remains the same.

Definition 3. The active set 𝒜(𝐮𝑠(0)) of any feasible 𝐮𝑠(0) satisfying (20) is the subset of {1,2,…,𝑁𝑚} such that 𝑗∈𝒜(𝐮𝑠(0)) if and only if 𝐚T𝑗𝐮𝑠(0)=𝐛𝑗.

Next, we present a one-step active set algorithm to solve the quadratic program (19)-(20) partially.

Algorithm 4. Given an initial feasible 𝐮𝑠(0) and a working set 𝒲0⊂𝒜(𝐮𝑠(0)). Let the set of working constraints 𝐚T𝑗𝐮𝑠(0)≤𝐛𝑗,𝑗∈𝒲0,(21) be represented by 𝐀0𝐮𝑠(0)≤𝐛0,(22) where the inequality sign applies componentwise, that is, each row of 𝐀0, 𝐛0 represents a working constraint in (21). (1)Compute the gradient 𝐠0=𝐇𝐮𝑠(0)+𝐅𝑣(0),(23) and the null space of 𝐀0, denoted 𝐙0 by 𝐀0𝐙0=0.(24) If 𝐙T0𝐠0≈0, go to step  (5). (2)Compute a search direction 𝐰0=𝐙0𝐰0, where 𝐙T0𝐇𝐙0𝐰0+𝐙T0𝐠0=0.(25)(3)Let 𝛼0=min𝐮𝑗∈𝒜𝑠(0)⧵𝒲0s.t.𝐚T𝑗𝐰0>0.𝐛𝑗−𝐚T𝑗𝐮𝑠(0)𝐚T𝑗𝐰0.(26)(4)If 𝛼0≥1, go to step  (5). Otherwise, update 𝐮𝑠(0) to 𝐮∗𝑠(0) by 𝐮∗𝑠(0)=𝐮𝑠(0)+𝛼0𝐰0,(27) 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)Update 𝐮𝑠(0) to 𝐮∗𝑠(0) by 𝐮∗𝑠(0)=𝐮𝑠(0)+𝐰0.(28) Compute the Lagrange multiplier 𝜆0 from 𝐀T0𝜆0+𝐠0=0(29) 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.

Algorithm 4 can be interpreted as follows. It solves the equality-constrained quadratic program:min𝐮′𝑠(0)ğ½î…žğ‘ ,ğ½î…žğ‘ î‚¸ğ®âˆ¶=î…žğ‘ î‚¹(0)𝑣(0)T𝐅𝐇𝐅Tğ†ğ®î‚¹î‚¸î…žğ‘ î‚¹(0)𝑣(0)subjectto𝐚Tğ‘—ğ®î…žğ‘ (0)=𝐛𝑗,𝑗∈𝒲0,,(30) and then searches along the direction𝐰0=ğ®î…žğ‘ (0)−𝐮𝑠(0),(31) until it is either blocked by a constraint not in 𝒲0 (step  (4)) or the optimal ğ®î…žğ‘ (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 𝐮∗𝑠(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 𝐽𝑠 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 𝑇𝑠. Let ğ‘ğ‘Žğ‘‡ğ‘  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 𝐮∗𝑠(0), 𝒲∗(0) to 𝐮𝑠(ğ‘ğ‘Ž), 𝒲(ğ‘ğ‘Ž) by rotating the components of 𝐮∗𝑠(0) by an amount of ğ‘ğ‘Ž since it is supposed to be periodic:âŽ¡âŽ¢âŽ¢âŽ¢âŽ¢âŽ¢âŽ¢âŽ¢âŽ¢âŽ£ğ‘¢ğ‘ î€·ğ‘ğ‘Žî€¸ğ‘¢ğ‘ î€·ğ‘ğ‘Žî€¸â‹®ğ‘¢+1𝑠𝑁𝑝𝑢−1𝑠𝑁𝑝𝑢𝑠𝑁𝑝⋮𝑢+1ğ‘ î€·ğ‘ğ‘Ž+ğ‘ğ‘î€¸âŽ¤âŽ¥âŽ¥âŽ¥âŽ¥âŽ¥âŽ¥âŽ¥âŽ¥âŽ¦=âŽ¡âŽ¢âŽ¢âŽ¢âŽ¢âŽ¢âŽ¢âŽ¢âŽ£ğ‘¢âˆ’1âˆ—ğ‘ î€·ğ‘ğ‘Žî€¸ğ‘¢âˆ—ğ‘ î€·ğ‘ğ‘Žî€¸â‹®ğ‘¢+1∗𝑠𝑁𝑝𝑢−1∗𝑠𝑢(0)∗𝑠⋮𝑢(1)âˆ—ğ‘ î€·ğ‘ğ‘Žî€¸âŽ¤âŽ¥âŽ¥âŽ¥âŽ¥âŽ¥âŽ¥âŽ¥âŽ¦âˆ’1.(32) Obviously, Algorithm 4 will then continue to solve an equivalent quadratic program as long as 𝑣 strictly follows the exogenous dynamics (14). Hence we have the following convergence result.

Proposition 5. Algorithm 4 together with the rotation of the components of 𝐮∗𝑠(0) in (32) will solve the quadratic program (19)-(20) in finite time as long as 𝑣(𝑘) satisfies the exogenous dynamics (14).

Proof. From the argument above it is easy to see that as long as 𝑣(𝑘) 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 in 𝑣. Since constraint (20) is independent of 𝑣, an abrupt change in 𝑣 will not destroy the feasibility of the previously determined 𝐮𝑠 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 Algorithm 4. One is to update the working set 𝒲 followed by 𝐮𝑠, and the other is to update 𝐮𝑠 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 𝐮𝑠 on an “outdated” working set if 𝑣 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 𝑣 has changed, and so has the objective function (19), the original stationary 𝐮∗𝑠 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 𝑣 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 𝐮𝑠(0) to ğ®î…žğ‘ (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 Algorithm 4 takes time ğ‘ğ‘Žğ‘‡ğ‘  to carry out, then 𝒲∗0 updated in the current iteration will be used to optimize 𝐮𝑠(ğ‘ğ‘Ž) in the next iteration, ğ®ğ‘ î€·ğ‘ğ‘Žî€¸âŽ¡âŽ¢âŽ¢âŽ¢âŽ£ğ‘¢âˆ¶=ğ‘ î€·ğ‘ğ‘Žî€¸ğ‘¢ğ‘ î€·ğ‘ğ‘Žî€¸â‹®ğ‘¢+1ğ‘ î€·ğ‘ğ‘Ž+ğ‘ğ‘î€¸âŽ¤âŽ¥âŽ¥âŽ¥âŽ¦âˆ’1,(33) but the outcome of that optimization is ready only at 𝑘=2ğ‘ğ‘Ž, based on which the transient control 𝑢𝑡 is computed. Suppose that the transient subproblem takes one sampling interval to solve, then the transient subproblem at 𝑘=2ğ‘ğ‘Ž will update 𝑢(2ğ‘ğ‘Ž+1)=𝑢𝑠(2ğ‘ğ‘Ž+1)+𝑢𝑡(2ğ‘ğ‘Ž+1) (see Section 5 below for a more detailed discussion). Hence, the “time priority” of 𝑢𝑠 is 2ğ‘ğ‘Ž+1,2ğ‘ğ‘Ž+2,…,𝑁𝑝−1,0,1,…,2ğ‘ğ‘Ž and from this argument, we choose to drop the constraint in ℒ that is associated with the first 𝑢𝑠 in this sequence or to add the constraint in ℬ that is associated with the last 𝑢𝑠 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 in 𝑣(𝑘). 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 𝐮𝑠(𝑘),𝒲(𝑘).

4. Steady-State Subproblem: A Dynamic Policy Approach

The approach proposed in Section 3 optimizes 𝐮𝑠. Consequently, the number of (scalar) variables being optimized is proportional to 𝑁𝑝. 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 𝑁𝑝, 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 𝑤 and 𝑟 in Section 3. Specifically, we assume 𝑢𝑠 is also generated from a state-space model: ̂𝑣+̂𝑣̂𝑢=𝐴𝑣,(34)𝑠̂𝑣̂=𝐶𝑣,(35) where𝐴𝑁𝑝̂𝑣=𝐼.(36) This state-space model is designed a priori but the initial state ̂𝑣(0) will be optimized online. Obviously, the quadratic program (19)-(20) becomeŝ𝑣min(0)𝐽𝑠,𝐽𝑠̂∶=𝑣(0)𝑣(0)T𝐇𝐅𝐅T𝐆̂,̂𝐚𝑣(0)𝑣(0)(37)subjecttoT𝑗̂𝑣(0)≤𝐛𝑗,𝑗=1,2,…,𝑁𝑚,(38) where𝒪𝐇=T𝐇̂𝐚𝒪,𝐅=𝐅𝒪,𝑗=𝐚𝑗𝒪,𝑗=1,2,…,𝑁𝑚,(39) and 𝒪 is the observability matrixîâŽ¡âŽ¢âŽ¢âŽ¢âŽ£ğ¶Ì‚ğ‘£ğ¶Ì‚ğ‘£ğ´Ì‚ğ‘£â‹®ğ¶Ì‚ğ‘£ğ´ğ’ªâˆ¶=𝑁𝑝−1Ì‚ğ‘£âŽ¤âŽ¥âŽ¥âŽ¥âŽ¦.(40) The number of variables in this new quadratic program is the dimension of ̂𝑣(0), denoted by 𝑛̂𝑣. If 𝐴̂𝑣 is constructed from the method of Fourier decomposition described in Section 3, Shannon’s sampling theorem implies that a sufficiently large but finite 𝑛̂𝑣 will guarantee a full reconstruction of the original optimization variable 𝐮𝑠(0). On the other hand, a smaller 𝑛̂𝑣 restricts the search to a lower dimensional subspace of 𝐮𝑠(0) and hence the optimization is easier but suboptimal.

One natural choice of the dynamics 𝐴̂𝑣 is to make 𝐴̂𝑣=𝐴𝑣 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 𝑢𝑠 may exceed that of 𝑤 and 𝑟 in order to achieve the true optimal performance. However, we could have over-designed the exogenous model (14)–(17) to include more harmonics in 𝐴𝑣 than necessary for 𝑤 and 𝑟, thus making the choice 𝐴̂𝑣=𝐴𝑣 here not so conservative. The simulation results in Section 6 will demonstrate both cases.

It remains to choose the matrix 𝐶̂𝑣 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., 𝑦𝑠=𝑟 or 𝑒𝑠=0) can be achieved by solving the regulator equation:𝑋𝐴𝑣=𝐴𝑋+𝐵1𝐶𝑤+𝐵2𝐶𝑈,𝑟=𝐶𝑋+𝐷1𝐶𝑤+𝐷2𝑈,(41) for the matrices 𝑋, 𝑈, and then let𝑢𝑠=𝑈𝑣,(42) which also implies𝑥𝑠=𝑋𝑣.(43) Therefore, to recover the optimal (or perfect) solution in the linear case when 𝑢𝑠 does not hit any constraint, the state-space model of 𝑢𝑠 may be chosen as̃𝑣+=𝐴𝑣̃𝑢𝑣,𝑠̃=𝑈𝑣.(44) However, this state-space model is not guaranteed to be observable. When it is not, the resulting 𝐇 in (37) will become semidefinite instead of strictly positive definite. To overcome this, we suggest to perform an orthogonal state transformation∗̂𝑣̃=𝑇𝑣,(45) to bring (44) to the Kalman decomposition form∗̂𝑣+=̂𝑣∗̂𝑣,𝑢∗∗0𝐴𝑠=̂𝑣∗̂𝑣,0𝐶(46) and hence obtain a reduced-order model to become (34)-(35). It is easy to verify that 𝐴𝑁𝑝̂𝑣=𝐼 since ̂𝑣∗∗0𝐴(47) is upper block-triangular and ̂𝑣∗∗0𝐴𝑁𝑝=𝑇𝐴𝑣𝑇−1𝑁𝑝=𝑇𝐴𝑁𝑝𝑣𝑇−1=𝐼.(48)

Certainly, the discussion above only provides a suggestion of how to choose the state-space model for 𝑢𝑠, 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 of ̂𝑣(0) rather than 𝐮𝑠(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 𝐮∗𝑠(0) to obtain 𝐮𝑠(ğ‘ğ‘Ž), we obtain ̂𝑣(ğ‘ğ‘Ž) by ğ´ğ‘ğ‘ŽÌ‚ğ‘£Ì‚ğ‘£âˆ—(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 Section 2. 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 𝐽𝑡 in (11), in actual implementation we shall optimize 𝑢𝑡(1),𝑢𝑡(2),…,𝑢𝑡(𝑁𝑢) based on 𝑥𝑡(1) at time 𝑘=0, instead of 𝑢𝑡(0),𝑢𝑡(1),…,𝑢𝑡(𝑁𝑢−1) based on 𝑥𝑡(0). The unknown 𝑥𝑡(1) can be projected from the known variables and system dynamics. After the optimization is carried out to obtain 𝑢𝑡(1),𝑢𝑡(2),…,𝑢𝑡(𝑁𝑢), the control input to be executed at 𝑘=1 will be 𝑢(1)=𝑢𝑠(1)+𝑢𝑡(1). Bear in mind that all these calculations can only be based on the best knowledge of the signals at 𝑘=0.

Figure 1 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 𝑢(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, 𝑣(0) is provided by the observer described in Section 3, 𝑥(0) is a measurement of the current plant state, and 𝑢(0) is the control input calculated by the algorithm at 𝑘=−1. Note that to compute 𝑢𝑡(1), the values of 𝑢𝑠(1),𝑢𝑠(2),…,𝑢𝑠(𝑁𝑢) are needed to form the constraints for the transient quadratic program since the original linear inequality constraints apply to 𝑢(𝑘)=𝑢𝑠(𝑘)+𝑢𝑡(𝑘). On the other hand, 𝑥𝑠(1) can be written as a linear function of 𝑣(1) and 𝐮𝑠(1) (or ̂𝑣(1)) explicitly. Finally, the steady-state subproblem requires a computational time of ğ‘ğ‘Žğ‘‡ğ‘ , implying that the solution 𝐮𝑠(0) provided by the steady-state subproblem at 𝑘=0 is based on a measurement of 𝑣(𝑘) at some 𝑘 between −2ğ‘ğ‘Ž+1 and âˆ’ğ‘ğ‘Ž. So in the worst case, 𝑢(1) is based on some information as old as 𝑣(−2ğ‘ğ‘Ž+1), which corresponds to a worst-case delay of 2ğ‘ğ‘Žğ‘‡ğ‘ . For instance, if ğ‘ğ‘Ž=1, the control 𝑢(𝑘) is computed from a measurement of 𝑥(𝑘−1), 𝑣(𝑘−1), and 𝑣(𝑘−2).

Remark 10. Although we said in Section 2 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 𝑢𝑡(𝑘) may have been totally consumed by 𝑢𝑠(𝑡) when the latter hits a constraint. For example, the simulations to be discussed in Section 6 have the constraint ||||=||𝑢𝑢(𝑘)𝑠(𝑘)+𝑢𝑡||(𝑘)≤1.(49) If 𝑢𝑠(𝑘) already saturates at ±1, one side of 𝑢𝑡(𝑘) is lost but that could be the only side to make the reduction of 𝐽𝑡 possible, thus forcing 𝑢𝑡(𝑘) to zero. So the input constraint does not only restrict the magnitude, but also the time instant to apply the transient control 𝑢𝑡. 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: ̇𝑥=𝐴𝑥+𝐵1𝑤+𝐵2𝑢,𝑦=𝐶𝑥+𝐷1𝑤+𝐷2âŽ¡âŽ¢âŽ¢âŽ¢âŽ¢âŽ¢âŽ¢âŽ¢âŽ¢âŽ¢âŽ¢âŽ£ğ‘¢,𝐴=−𝑅𝑙𝐿𝑙00−1𝐿𝑙−1𝐿𝑙0−𝑅se𝐿se0−1𝐿se000−𝑅sh𝐿sh0−1𝐿sh1𝐶se1𝐶se1000𝐶sh01𝐶sh⎤⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦,𝐵001=⎡⎢⎢⎢⎢⎢⎢⎣1𝐿𝑙00000000−1𝐶sh⎤⎥⎥⎥⎥⎥⎥⎦,𝐵2=⎡⎢⎢⎢⎢⎢⎢⎣V00dc2𝐿se00Vdc2𝐿sh⎤⎥⎥⎥⎥⎥⎥⎦,,𝐷0000𝐶=00001100001=0000,𝐷2=.0000(50) The exogenous input 𝑤 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 𝑦 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 𝑢 is composed of two switching signals across the voltage source inverters (VSIs), both of which are required to satisfy the bounds −1≤𝑢≤1. The general control objective is to maintain 𝑦 to the desired waveforms despite possible fluctuations in 𝑤 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 𝑇𝑠=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 𝑢𝑡 depends on 𝑢𝑠, the transient control 𝑢𝑡(𝑘+1) is computed from 𝐮𝑠(𝑘) 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 Sections 3 and 4; 𝑁𝑝 may be chosen to represent half of the period instead of the whole period, satisfying𝑥𝑠(0)=−𝑥𝑠𝑁𝑝,(51) instead of (6), and𝐴𝑁𝑝𝑣𝐴=−𝐼,𝑁𝑝̂𝑣=−𝐼,(52) instead of (17) and (36), with 𝜔𝑇𝑠×𝑁𝑝=𝜋 instead of 2𝜋. The rotation operation in (32) should also becomeâŽ¡âŽ¢âŽ¢âŽ¢âŽ¢âŽ¢âŽ¢âŽ¢âŽ¢âŽ£ğ‘¢ğ‘ î€·ğ‘ğ‘Žî€¸ğ‘¢ğ‘ î€·ğ‘ğ‘Žî€¸â‹®ğ‘¢+1𝑠𝑁𝑝𝑢−1𝑠𝑁𝑝𝑢𝑠𝑁𝑝⋮𝑢+1ğ‘ î€·ğ‘ğ‘Ž+ğ‘ğ‘î€¸âŽ¤âŽ¥âŽ¥âŽ¥âŽ¥âŽ¥âŽ¥âŽ¥âŽ¥âŽ¦=âŽ¡âŽ¢âŽ¢âŽ¢âŽ¢âŽ¢âŽ¢âŽ¢âŽ£ğ‘¢âˆ’1âˆ—ğ‘ î€·ğ‘ğ‘Žî€¸ğ‘¢âˆ—ğ‘ î€·ğ‘ğ‘Žî€¸â‹®ğ‘¢+1∗𝑠𝑁𝑝−1−𝑢∗𝑠(0)−𝑢∗𝑠⋮(1)âˆ’ğ‘¢âˆ—ğ‘ î€·ğ‘ğ‘Žî€¸âŽ¤âŽ¥âŽ¥âŽ¥âŽ¥âŽ¥âŽ¥âŽ¥âŽ¦âˆ’1.(53) 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 Table 1 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 𝑢 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.

Component Vdc 𝐿 s e 𝐿 s h 𝐶 s e 𝐶 s h

Value320 V 5.0 mH 1.2 mH 10  𝜇 F 20  𝜇 F

Component 𝑅 𝑙 𝐿 𝑙 𝑅 s e 𝑅 s h

Value 0.01  Ω 1.0 mH 0.01  Ω 0.01  Ω

The simulation scenario is summarized in Figure 2. 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 𝑤-dependent. The waveforms of both 𝑤 and 𝑟 are shown in Figure 2. The simulation scenario is designed such that the steady-state control 𝑢𝑠 is not saturated at the beginning. At 𝑡=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 𝑢𝑠 into a slightly saturated situation. At 𝑡=1.0s, the load demand increases, causing the reference of the supply current to increase again, and 𝑢𝑠 will become deeply saturated. At 𝑡=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 𝑢𝑠 required will just be within the bounds of ±1, thus leaving the saturation region to return to the linear region. So, in short, 𝑢𝑠 is expected to experience “linear → slightly saturated → deeply saturated →linear but nearly saturated” in this simulation scenario.

To evaluate the performance of our algorithms, we compare them to two other cases. In the first case, instead of Algorithm 4, the complete quadratic program is solved in each iteration of the steady-state subproblem every ğ‘ğ‘Ž-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 𝐮𝑠(0) in the steady-state subproblem is supposed to be𝐮𝑠(âŽ¡âŽ¢âŽ¢âŽ¢âŽ£ğ‘ˆ0)=𝑈𝐴𝑣⋮𝑈𝐴𝑁𝑝𝑣−1âŽ¤âŽ¥âŽ¥âŽ¥âŽ¦ğ‘£(0),(54) where 𝑈 is the solution to the regulator equation (41). The transient subproblem can also be solved by𝑢𝑡=𝐹𝑥𝑡,(55) where 𝐹 is the optimal state-feedback gain minimizing the transient quadratic cost âˆžî“ğ‘˜=0𝑦𝑡(𝑘)T𝑄𝑦𝑡(𝑘)+𝑢𝑡(𝑘)T𝑅𝑢𝑡(𝑘).(56) However, the combined input 𝑢 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 𝐽𝑠 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 Section 5 will be in force, where 𝑢(𝑘+1) instead of 𝑢(𝑘) will be optimized and the steady-state control 𝐮𝑠 is only updated every ğ‘ğ‘Ž-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.

Figure 3 plots the steady-state cost 𝐽𝑠 of our first approach in Section 3 based on receding horizon quadratic programming together with the costs in the other two cases. ğ‘ğ‘Ž is assumed to be 3 in this simulation. The transient subproblem has a control horizon of 𝑁𝑢=5, corresponding to 10 scalar variables and 20 scalar inequality constraints. On the other hand, the steady-state subproblem has 𝑁𝑝=50 corresponding to 100 variables and 200 constraints. As shown in Figure 3, all 𝐽𝑠 are zero prior to 𝑡=0.5s and should also settle down to zero after 𝑡=1.5s. It is observed that the transient response of our 𝐽𝑠 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 𝐽𝑠 value when a saturation occurs, due to its ignorance of the control constraints. The exact values of 𝐽𝑠 just prior to 𝑡=1.0s and 𝑡=1.5s are summarized in Table 3, which are about 6-7 times the 𝐽𝑠 values of our algorithm.

Steady-state cost 𝐽 𝑠 𝑁 ğ‘Ž Number of
at 𝑡 < 1 . 0 s at 𝑡 < 1 . 5 sVariablesConstraints

MVR 0.1156 53.276
MPC, receding horizon quadratic programming 0.0178 8.7527 3100 ( 𝑁 𝑝 × 2 )200 ( 𝑁 𝑝 × 2 × 2 )
MPC, dynamic policy (up to 9th harmonics) 0.0225 9.0273 1 𝑛 ̂ 𝑣 = 2 0 ( 5 × 2 × 2 )200 ( 𝑁 𝑝 × 2 × 2 )
MPC, dynamic policy (up to 29th harmonics) 0.0182 8.7640 2 𝑛 ̂ 𝑣 = 6 0 ( 1 5 × 2 × 2 )200 ( 𝑁 𝑝 × 2 × 2 )
Transient subproblem10 ( 𝑁 𝑢 × 2 )20 ( 𝑁 𝑢 × 2 × 2 )

Figure 4 zooms into the first component of the control input 𝑢. The top two plots draw attentions to the steady-state control 𝑢𝑠 (𝑢 is essentially 𝑢𝑠 just prior to 𝑡=1.0s and 𝑡=1.5s). It is observed that the differences in 𝑢𝑠 between MVR and MPC are very subtle. Compared with the MVR, the MPC 𝑢𝑠 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 𝐽𝑠 down by almost one order of magnitude. This demonstrates how nontrivial the optimal 𝑢𝑠 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 (𝑡>1.0s and 𝑡>1.5s), the transient control 𝑢𝑡 is “disabled” due to the saturation of the steady-state control 𝑢𝑠. Note that we are labeling the dashed blue curve as “𝑢𝑠 used to compute 𝑢𝑡” since it is slightly different from the actual 𝑢𝑠. For instance, 𝑢𝑡(𝑘+1) is computed from the knowledge of 𝑢𝑠 at time 𝑘, which is not exactly the same as 𝑢𝑠(𝑘+1). Obviously, 𝑢𝑡 is not just disabled whenever 𝑢𝑠 saturates. It happens only when the desired direction of 𝑢𝑡 violates the active constraint.

Next, let us look at the performance of our second MPC approach in Section 4 based on dynamic policy. Odd harmonics up to the 9th order are included in 𝐴̂𝑣 resulting in a total of 𝑛̂𝑣=20 variables. See Table 3. Since the number of variables is much lower than that of the first approach, we assume ğ‘ğ‘Ž=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 𝐽𝑠 is plotted in Figure 5. Note that the MVR curve exhibits a slightly different transient from Figure 3 since their ğ‘ğ‘Ž 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 ğ‘ğ‘Ž but also a smaller-sized quadratic program overall. However, the drawback is a slightly suboptimal 𝐽𝑠, as indicated in Table 3.

As mentioned in Section 4, it is possible to over-design 𝐴𝑣, and hence 𝐴̂𝑣, so that the optimal 𝐽𝑠 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 𝑤, we may include odd harmonics up to the 29th order in 𝐴𝑣 and 𝐴̂𝑣. The results are also recorded in Table 3, and we see that this 𝐽𝑠 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.


This work is supported by HKU CRCG 201008159001.


  1. J. M. Maciejowski, Predictive Control with Constraints, Prentice-Hall, New Jersey, NY, USA, 2002.
  2. S. J. Qin and T. A. Badgwell, “A survey of industrial model predictive control technology,” Control Engineering Practice, vol. 11, no. 7, pp. 733–764, 2003. View at: Publisher Site | Google Scholar
  3. P. Whittle, Optimization Over Time, Wiley, New York, NY, USA, 1982.
  4. W. H. Kwon and S. Han, Receding Horizon Control, Springer, New York, NY, USA, 2005.
  5. A. Bemporad, M. Morari, V. Dua, and E. N. Pistikopoulos, “The explicit linear quadratic regulator for constrained systems,” Automatica, vol. 38, no. 1, pp. 3–20, 2002. View at: Publisher Site | Google Scholar
  6. A. Bemporad, F. Borrelli, and M. Morari, “Min-max control of constrained uncertain discrete-time linear systems,” IEEE Transactions on Automatic Control, vol. 48, no. 9, pp. 1600–1606, 2003. View at: Publisher Site | Google Scholar
  7. A. Bemporad and C. Filippi, “Suboptimal explicit receding horizon control via approximate multiparametric quadratic programming,” Journal of Optimization Theory and Applications, vol. 117, no. 1, pp. 9–38, 2003. View at: Publisher Site | Google Scholar
  8. Z. Wan and M. V. Kothare, “Robust output feedback model predictive control using off-line linear matrix inequalities,” Journal of Process Control, vol. 12, no. 7, pp. 763–774, 2002. View at: Publisher Site | Google Scholar
  9. Z. Wan and M. V. Kothare, “An efficient off-line formulation of robust model predictive control using linear matrix inequalities,” Automatica, vol. 39, no. 5, pp. 837–846, 2003. View at: Publisher Site | Google Scholar
  10. M. Åkerblad and A. Hansson, “Efficient solution of second order cone program for model predictive control,” International Journal of Control, vol. 77, no. 1, pp. 55–77, 2004. View at: Publisher Site | Google Scholar
  11. C. V. Rao, S. J. Wright, and J. B. Rawlings, “Application of interior-point methods to model predictive control,” Journal of Optimization Theory and Applications, vol. 99, no. 3, pp. 723–757, 1998. View at: Publisher Site | Google Scholar
  12. J. Mattingley, Y. Wang, and S. Boyd, “Receding horizon control: automatic generation of high-speed solvers,” IEEE Control Systems Magazine, vol. 31, no. 3, pp. 52–65, 2011. View at: Google Scholar
  13. M. Cannon and B. Kouvaritakis, “Optimizing prediction dynamics for robust MPC,” IEEE Transactions on Automatic Control, vol. 50, no. 11, pp. 1892–1897, 2005. View at: Publisher Site | Google Scholar
  14. A. Richards, K.-V. Ling, and J. Maciejowski, “Robust multiplexed model predictive control,” in Proceedings of the European Control Conference, Kos, Greece, July 2007. View at: Google Scholar
  15. K. V. Ling, W. K. Ho, Y. Feng, and B. Wu, “Integral-square-error performance of multiplexed model predictive control,” IEEE Transactions on Industrial Informatics, vol. 7, no. 2, pp. 196–203, 2011. View at: Google Scholar
  16. D. Li and Y. Xi, “Aggregation based robust model predictive controller for systems with bounded disturbance,” in Proceedings of the 7th World Congress on Intelligent Control and Automation, (WCICA '08), pp. 3374–3379, Chongqing, China, June 2008. View at: Publisher Site | Google Scholar
  17. D. Li and Y. Xi, “Aggregation based closed-loop MPC with guaranteed performance,” in Proceedings of the 48th IEEE Conference on Decision and Control, (CDC/CCC '09), pp. 7400–7405, Shanghai, China, December 2009. View at: Publisher Site | Google Scholar
  18. J. B. Rawlings and K. R. Muske, “The stability of constrained receding horizon control,” IEEE Transactions on Automatic Control, vol. 38, no. 10, pp. 1512–1516, 1993. View at: Publisher Site | Google Scholar
  19. K. H. Kwan, Y. C. Chu, and P. L. So, “Model-based H∞ control of a unified power quality conditioner,” IEEE Transactions on Industrial Electronics, vol. 56, no. 7, pp. 2493–2504, 2009. View at: Publisher Site | Google Scholar
  20. S. G. Nash and A. Sofer, Linear and Nonlinear Programming, McGraw-Hill, New York, NY, USA, 1996.
  21. J. Nocedal and S. J. Wright, Numerical Optimization, Springer, New York, NY, USA, 2nd edition, 2006.
  22. A. Gautam, Y.-C. Chu, and Y. C. Soh, “Robust MPC with optimized controller dynamics for linear polytopic systems with bounded additive disturbances,” in Proceedings of the 7th Asian Control Conference, (ASCC '09), pp. 1322–1327, Hong Kong, China, August 2009. View at: Google Scholar
  23. A. Gautam, Y.-C. Chu, and Y. C. Soh, “Optimized dynamic policy for receding horizon control of linear time-varying systems with bounded disturbances,” IEEE Transactions on Automatic Control. In press. View at: Google Scholar
  24. E. J. Davison, “The robust control of a servomechanism problem for linear time-invariant multivariable systems,” IEEE Transactions on Automatic Control, vol. 21, no. 1, pp. 25–34, 1976. View at: Google Scholar
  25. B. A. Francis and W. M. Wonham, “The internal model principle of control theory,” Automatica, vol. 12, no. 5, pp. 457–465, 1976. View at: Google Scholar
  26. B. A. Francis, “The linear multivariable regulator problem,” SIAM Journal on Control Optimization, vol. 15, no. 3, pp. 486–505, 1977. View at: Google Scholar

Copyright © 2012 Yun-Chung Chu and Michael Z. Q. Chen. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

More related articles

 PDF Download Citation Citation
 Download other formatsMore
 Order printed copiesOrder

Related articles