Constrained Model Predictive Control using Feedforward Steady-State
Target Optimization tracking approach for an Isolated Power System
Paulo R. Loma Marconi
1
Abstract A constrained Linear Quadratic Model Predictive
Model Control (LQ-MPC) with offset-free tracking and dis-
turbance rejection using the Feedforward Steady-State Target
Optimization (SSTO) approach is presented in this paper.
The MPC control law uses Finite-Receding-Horizon method
in dual-mode. Moreover, the stability and recursive feasibility
are guaranteed by off-line dual-mode and on-line optimization.
The results are presented by simulation.
I. STATE-SPACE MODEL OF THE SYSTEM
The isolated power system [1] to be studied is a linear
Single-Input Single-Output (SISO) model with the addition
of the disturbance in the output, Fig. 1.
Fig. 1: Isolated power system
The linear continuous-time state-space model is as follows,
˙ω
˙p
m
˙p
v
=
D/M 1/M 0
0 1/T
t
1/T
t
1/(RT
g
) 0 1/T
g
ω
p
p
v
+
0
0
1/T
g
p
ref
+
1/M
0
0
p
d
f =
50 0 0
ω
p
p
v
+ p
d
(1)
where f is the output values, p
ref
is the input, and p
d
is the disturbance, the rest of the parameters can be consulted
at [1].
Thus, the Discrete Linear Time-Invariant (LTI) state-space
system is
x(k + 1) = A x(k) + B u(k) + B
d
d(k)
y(k) = C x(k) + D
d
d(k), k = 0, 1, 2, . . .
(2)
where x R
n
, u R
m
, y R
p
, and d R
q
are the states,
input, output, and disturbance vectors respectively. A, B, C,
are the state, input, and output matrices, and B
d
, D
d
are the
disturbance matrices in the input, and in the output.
1
The author is with the Department of Automatic Con-
trol Systems Engineering, The University of Sheffield, UK
prlomamarconi1@sheffiel.ac.uk
II. FORMULATION OF THE LQ-MPC TRACKING
PROBLEM
A. Feedforward tracking with SSTO
The Feedfoward tracking approach by SSTO requires a
target optimizer that calculates the steady-state values for
the states x
ss
and the input u
ss
given a reference r in the
presence of a disturbance d, Fig. 2.
Fig. 2: MPC and SSTO model
B. Cost function
Because it is desired that x x
ss
and u u
ss
, new
variables can be defined as ξ x x
ss
and v u u
ss
,
where ξ 0 and v 0.
Therefore, the constrained LQ-MPC problem is to mini-
mize the cost function
IP
N
(v(k)) : min
v (k)
N1
X
j=0
(ξ
|
(k + j|k) Q ξ(k + j|k)
+ v
|
(k + j|k) R v(k + j|k))
+ ξ
|
(k + N|k) P ξ(k + N|k)
(3)
subject to, for j = 0, 1, 2, . . . , N 1
ξ(k + 1 + j|k) = A ξ(k + j|k) + B v(k + j|k)
ξ(k|k) = ξ(k)
P
x
ξ(k + j|k) q
x
P
x
x
ss
(4)
P
u
v(k + j|k) q
u
P
u
u
ss
(5)
P
x
N
ξ(k + N|k) ˜q
x
N
(6)
where Q and R are weight matrices, P is the terminal cost
matrix, (4), (5), (6) are the linear inequality (polyhedra)
constraints, and N is the finite-horizon value. Note that
B
d
, D
d
and d are not present because the SSTO will manage
the disturbance rejection.
C. Defining the inequality constraints
Given |p
ref
| = |u| u
min,max
, and |f| = |y|
y
min,max
, for y = C x,
C
C
x(k + j|k)
y
max
y
min
I
I
u(k + j|k)
u
max
u
min
can be written in terms of ξ and v as follows,
C
C
| {z }
P
x
ξ(k + j|k)
y
max
y
min
| {z }
q
x
C
C
| {z }
P
x
x
ss
I
I
|{z}
P
u
v(k + j|k)
u
max
u
min
| {z }
q
u
I
I
|{z}
P
u
u
ss
(7)
and using a proper deadbeat mode-2 terminal inequality
constraint to ensure a fast convergence to zero,
I
n×n
P
x
P
u
K
(A + BK
)
0
.
.
.
(A + BK
)
N1
| {z }
P
x
N
ξ(k + N|k)
1
n
q
x
P
x
x
ss
q
u
P
u
u
ss
| {z }
˜q
x
N
(8)
where K
is the deadbeat mode-2 gain calculated using
Linear Quadratic Regulator (LQR) method.
D. Target equilibrium optimization under constraints (SSTO)
The offset-free equilibrium points (x
ss
, u
ss
) that satisfies
the constraints exits if only if x
ss
= Ax
ss
+ Bu
ss
+ B
d
d
and y
ss
= Cx
ss
+ D
d
d = r.
Thus, the SSTO problem is
min
{x
ss
,u
ss
}
kC x
ss
+ D
d
d rk + ρku
ss
k (9)
subject to,
I A B
C 0
| {z }
A
eq
=T
x
ss
u
ss
|{z}
x
eq
=
B
d
r D
d
d
| {z }
b
eq
P
x
0
0 P
u
| {z }
A
neq
x
ss
u
ss
|{z}
x
neq
q
x
q
u
|{z}
b
neq
(10)
III. SOLVING THE LQ-MPC OPTIMIZATION
PROBLEM
The solution of the LQ-MPC problem (3) in the compress
form can written as,
min
v (k)
1
2
v
|
(k) H v(k) + ξ
|
(k) L
|
v(k) + ξ
|
(k) M ξ(k)
(11)
subject to
P
c
v(k) q
c
+ S
c
ξ(k) (12)
where H, L, M are the cost matrices given Q and
R. P
c
, q
c
, S
c
are the stacked inequality constraints for
P
x
, P
u
, P
x
N
, q
u
, q
x
, ˜q
x
N
given the predictions matrices x =
F x(k) + Gu(k), [1]. Note that the prediction matrices are
constructed using the states and input without the tracking
model.
The solution of (11) is the optimal control input sequence,
v
= {v
(k|k), v
(k + 1|k), . . . , v
(k + N 1|k)}
which can be solved numerically.
Because MPC uses the receding-horizon principle, only
the first optimized control input is needed
v(k) = v
(k|k) = K
N
(ξ(k))
but the real control input applied to the system is
u(k) = u
ss
+ v
(k|k) (13)
A. Stability
In order to ensure stability it is necessary to compute a
terminal cost matrix P that satisfies the Lyapunov equation
(A + B K)
|
P (A + B K) P + (Q + K
|
R K) = 0 (14)
which solution can be computed numerically. A good ap-
proach is to use a deadbeat mode-2 control gain K that
converges the states to the origin x = 0 as long as (A+BK)
is stable, which means that the eigenvalues λ are inside the
unit circle |λ| < 1.
IV. ALGORITHM FOR THE IMPLEMENTATION
Now that the solution of the LQ-MPC is established and
the stability can be guaranteed, there are basic conditions to
meet such as,
(A, B) is stabilizable which can be proved by the
controllability, in Matlab rank(ctrb(A,B)).
Q 0, can be C
|
C
R 0, changed later for tuning.
The following steps reflect the Matlab implementation
algorithm
Off-line mode:
1. Discretization of the model (1) using
zero holder (zoh), in Matlab will be
sys=c2d(idss(Ac,Bc,C,D,Bd,x0,0),Ts),
where T
s
is the sampling time.
2. Construct the target equilibrium constraints (10)
3. Define the optimization equation as z=[xss;uss]
and fun=@(z) C(1)
*
z(1)+Dd
*
d-r+rho
*
z(4)
4. Solve the optimization problem for SSTO using
f=fmincon(fun,z0,Aneq,bneq,Aeq,beq) in
Matlab, where z0=[r,0,0,0], fun is (9), and
f=[xss,uss] is the optimized vector solution.
5. Construct the LQ-MPC inequality constraints (11), and
compute the deadbeat mode-2 terminal constraint using
K=-dlqr(A,B,Q,R).
6. Construct the prediction matrices F, G using
[F,G]=predict_mats(A,B,N) and the
constraint matrices (12) using [Pc,qc,Sc]=...
constraint_mats(F,G,Pu,qu,Px,qx,PxN,qxN);
7. Compute the deadbeat mode-2 control Kl using
Kl=-acker(A,B,[0,0,0])
8. Calculate the Lyapunov solution of (14) with
P=dlyap((A+B
*
K)’,Q+K’
*
R
*
K), which gives a
terminal cost matrix P that guaranties stability.
9. Finally before the on-line mode is initiated, set the
initial conditions as x=[0;0;0;]
On-line mode: iteration for k steps
1. Recall steps 3-4 of the off-line mode for the SSTO
optimization but at each iteration z0=[x;us(1,k)],
where us(1,k) is the first optimized control input.
2. Calculate ξ(k) and v(k) using xi(:,k)=x-xss
v(:,k)=us(:,k)-uss.
3. Compute the constraints (10) with the new x
ss
and u
ss
4. Compute (12) like it was done in the off-line mode
step 6.
5. Store the states xs(:,k)=x for the iteration process.
6. Solve the LQ-MPC optimization problem
(11) without using the matrix M as follows,
[v,fval,flag]=quadprog(H,L
*
xi(:,k),...
Pc,qc+Sc
*
xi(:,k)) where v is the optimal control
input v(k)
7. Check feasibility using flag, if it is < 1 the solution
is infeasible.
8. Use (13) as the real control input applied to the sys-
tem, us(:,k)=v(1)+uss using only the first
optimized value (receding principle).
9. Close the loop (2), x=A
*
x+B
*
us(1,k)+Bd
*
d
y=C
*
xs+Dd
*
d
10. Calculate the cost value of (11) using the matrix M .
11. Return to step 1.
V. DESIGN REQUIREMENTS
For the simulation the following specifications are de-
manded,
Reference r = 0.3
Frequency settles to |f | = |y| 0.01 within 2 [s] of
the disturbance initiaion.
|p
ref
| = |u| 0.5, and |f| = |y| 0.5.
Zero steady-state error in the output.
VI. SIMULATION AND RESULTS
The tuning and design parameters, as well as the simula-
tion value results are shown in TableI
Fig. 3a shows that the input constraint is satisfied, and in
Fig. 3b it can be seen that the zero offset-free tracking is also
satisfied rejecting disturbances at different times. The low
value of R is to minimize the overshoot as much as possible
which is 0.6%, the best Q is C
|
C where the settling time
t
s
is around 0.7. Also, the selection of N = 3 was lowest
possible to increase performance, Fig. 3c.
Parameter Value
Q C
|
C
R 0.08
N 3
d
1
0.5 at t = 30[s]
d
2
0.2 at t = 60[s]
r 0.3
x
0
[0; 0; 0]
|
t
s
0.7[s]
overshoot 0.6%
TABLE I: Tuning and design parameters
(a) u(k)
(b) y(k) vs r(k)
(c) y(k) vs r(k)
Fig. 3: Input and output
Fig. 4 shows that the constraints are satisfied for the three
states. And the cost value shown in Fig. 5 is relatively low
due the low value of R.
(a) States vs step (k)
(b) x1(k) vs x2(k)
(c) x1(k) vs x3(k)
(d) x2 (k) vs x3(k)
Fig. 4: States and phase plots
Fig. 5: Cost value J(k) function
One of the drawbacks of the SSTO approach is that it does
not have robustness because it needs a perfect model of the
plant. Also, when d > 0.5, the solution is infeasible, in order
to handle bigger disturbances it is necessary to change the
cost function, [2] shows an alternative solution.
REFERENCES
[1] P. Trodden, Lecture Notes ACS616 2018/19
[2] S. Dughman, J.A. Rossiter, A survey of guaranteeing feasibility and
stability in MPC during target changes 9th International Symposium
on Advanced Control of Chemical Processes, June 2015.