![](bg5.png)
where s ∈ R
n
+
is the slack variable and w ∈ R
n
+
∪ ∞ is
its weight (Giselsson, 2013; Kouzoupis, 2014). If the
QP in (4) has a feasible solution and if the weight w
is big enough, the solution of the new problem is the
solution of the original QP and s = 0. For smaller w
or for infeasible QP, s has non-zero components.
A way to efficiently implement soft constraints
with linear cost of constraint violation can be seen by
comparing (5) to (30). If w is taken to be the upper
bound for Lagrange multipliers v, the dual proximal
gradient method solves the soft-constrained problem.
The components of the dual residual corresponding to
the violated soft constraints are 0, thus violated soft
state constraints appear among inactive constraints
when calculating the local rate of decrease of resid-
ual.
5 EXAMPLE
A discrete-time form of the AFTI-16 benchmark
model as in (Giselsson, 2013) has the system matri-
ces
A=
0.9993 −3.0083 −0.1131 −1.6081
−0.0000 0.9862 0.0478 0.0000
0.0000 2.0833 1.0089 −0.0000
0.0000 0.0526 0.0498 1.0000
B=
−0.0804 −0.6347
−0.0291 −0.0143
−0.8679 −0.0917
−0.0216 −0.0022
in (1). The constraints are:
X =
x ∈ R
4
;−0.5 − s
1
≤ x
2
≤ 0.5 + s
1
,
−100 − s
2
≤ x
4
≤ 100 + s
2
}
U =
u ∈ R
2
;−25 ≤ u
1
≤ 25,
−25 ≤ u
2
≤ 25
}
. (32)
The constraints on the components of x are soft, the
linear weight on the components of the slack is w =
[10
5
, 10
5
]
T
, and w
T
s(k) is added to the sum term in
(2). The cost matrices are
Q = diag(10
−4
, 10
2
, 10
−3
, 10
2
),
R = diag(10
−2
, 10
−2
). (33)
The reference x
r
is 0 in all components except for the
first 50 time steps of the simulation, during which x
4
is 10. Initial state at the beginning of the simulation is
x(0) = [0, 0, 0, 0]
T
.
A family of QP in condensed form (4) correspond-
ing to the MPC problem is formed for N = 10 using
QPgen (QPgen, 2014; Giselsson and Boyd, 2014).
The matrices C and H are constant while the vectors
b and c are dependent on the system state and the ref-
erence. The QPs are 20-dimensional (10 x 2 compo-
nents of u) with 40 lines in C (limits on 10 x 2 input
signals and 10 x 2 state components).
The preconditioning diagonal matrix E is chosen
so as to minimize the condition number of the non-
singular part of M while setting the highest eigen-
value of M to 1. QPgen finds E = diag(10.4796,
3.5413, 9.9973, 10.0080, 9.9987, 10.0005, 10.0000,
10.0037, 9.9990, 9.9997, 10.0001, 10.0033, 9.9989,
9.9979, 10.0003, 10.0036, 9.9999, 9.9965, 9.9972,
10.0034, 0.2058, 0.0918, 0.1003, 0.1000, 0.1007,
0.1001, 0.1005, 0.1000, 0.1004, 0.1001, 0.1007,
0.1000, 0.1004, 0.1001, 0.1009, 0.0999, 0.1004,
0.1000, 0.1013, 0.1000). The model is initially simu-
lated in closed loop with 10
6
iterations in every time
step for 100 time steps. The system state x is recorded
at every time step and used as the input initial state
for observing convergence of the resulting QPs. Al-
gorithm behaviour is analysed in all time steps for 1
to 3000 iterations.
Figure 1 shows the active sets through iterations
for the first 10 time steps. White columns are delim-
iters of time steps.
In Figure 2, convergence through the active set
changes in time step 1 is shown graphically. Firstly
we analyse the behaviour with the final ω. The
quadratic norm of the primal residual is 0.22805 in
the 300th iteration and 8.74325×10
−12
in the 1800th
iteration, so the difference gets multiplied by 0.98414
in each iteration. This gives the estimate for the small-
est non-zero eigenvalue of the relevant M[ω|ω] of
1 − 0.98414 = 0.01586. In fact, 20 constraints are
active during the considered iterations (and no soft
ones are violated), 9 of them soft corresponding to
state constraints, and 11 hard corresponding to input
constraints. The matrix M[ω|ω] is non-singular and
its lowest eigenvalue is 0.01587, showing good agree-
ment with the numerical behaviour.
There is another longer interval in the first sam-
ple where the active set does not change between the
24th and 153rd iteration. The set has 22 elements, 9
of which are soft constraints and 13 are hard; no con-
straint is violated. If the dual residual is transformed
into the eigenspace of M[ω|ω], one expects the com-
ponents to decay in proportion with their correspond-
ing eigenvalues of M. These components are plotted
in Figure 3. It can be seen that the 22 lines decay
with different slopes in the relevant region and that the
ones listed higher in the legend decay slower. They
are listed in the figure in ascending order with respect
to eigenvalues (the first 2 correspond to the eigenval-
ues that are 0 and stay constant). Similarly, the active
set is constant from the 161st iteration on, has 20 ele-
ICINCO 2017 - 14th International Conference on Informatics in Control, Automation and Robotics
58