Note
Go to the end to download the full example code.
MPC controller for PWA systems#
This example demos how an MPC controller can be built for a simple system with piecewise
affine dynamics using csnlp.wrappers.PwaMpc. We assume some knowledge of MPC
is already present; otherwise, refer to other optimal control examples. For more details
on the subject of control for PWA and hybrid systems , see [2]
and [3].
Briefly, PWA systems are systems whose dynamics are described by a set of affine systems that are active in different regions of the state space. For region \(i\), the dynamics are described as
\[x_+ = A_i x + B_i u + c_i \quad \text{if} \quad S_i [x^\top, u^\top]^\top \leq T_i.\]
In the context of optimal control, with the proper procedure, these dynamics can be translated into a mixed-integer optimization problem. To do so, polytopic bounds on the state and input variables must be defined as
\[D [x^\top, u^\top]^\top \leq E.\]
The procedure to convert the PWA system into a mixed-integer optimization problem requires solving linear programs, whose number increases with the number of regions in the system. So, the fewer regions, the computationally lighter building the MPC is.
We start with the imports as usual. All MPC controller classes can be found in the
csnlp.wrappers module.
import casadi as cs
import matplotlib.pyplot as plt
import numpy as np
from csnlp import Nlp, wrappers
Setup#
Dynamics#
We consider a simple one-sided spring-mass-damper system with piecewise affine dynamics. In particular, the system has two modes: one where the spring is active and another where it is not, thus yielding two regions with different affine dynamics.
tau = 0.5 # sampling time for discretization
k1 = 10 # spring constant when one-sided spring active
k2 = 1 # spring constant when one-sided spring not active
damp = 4 # damping constant
mass = 10 # mass of the system
A_spring_1 = np.array([[1, tau], [-((tau * 2 * k1) / mass), 1 - (tau * damp) / mass]])
A_spring_2 = np.array([[1, tau], [-((tau * 2 * k2) / mass), 1 - (tau * damp) / mass]])
B_spring = np.array([[0], [tau / mass]])
From these matrices, we can define a sequence of csnlp.wrappers.PwaRegion
objects that represent the dynamics of the system in each region.
pwa_system = (
wrappers.PwaRegion(
# region dynamics
A=A_spring_1,
B=B_spring,
c=np.zeros(2),
# region domain
S=np.array([[1, 0, 0]]),
T=np.zeros(1),
),
wrappers.PwaRegion(
# region dynamics
A=A_spring_2,
B=B_spring,
c=np.zeros(2),
# region domain
S=np.array([[-1, 0, 0]]),
T=np.zeros(1),
),
)
Bounds#
In order for the PWA system to be converted to a mixed-integer optimization problem, we need to define the bounds of the system. In this case, we must impose polytopic bounds D @ [x; u] <= E on the states and the inputs as follows.
PWA MPC Controller#
The MPC controller is built in the same way as for other MPC classes. The main
difference is that now we must use the csnlp.wrappers.PwaMpc.set_pwa_dynamics
instead of csnlp.wrappers.Mpc.set_affine_dynamics to set the dynamics of the
system.
N = 10
mpc = wrappers.PwaMpc(nlp=Nlp[cs.SX](sym_type="SX"), prediction_horizon=N)
x, _ = mpc.state("x", 2)
u, _ = mpc.action("u")
mpc.set_pwa_dynamics(pwa_system, D, E)
mpc.constraint("state_constraints", D1 @ x - E1, "<=", 0)
mpc.constraint("input_constraints", D2 @ u - E2, "<=", 0)
mpc.minimize(cs.sumsqr(x) + cs.sumsqr(u))
mpc.init_solver({"record_time": True}, "bonmin") # "bonmin", "knitro", "gurobi"
Clp0006I 0 Obj 0 Dual inf 0.0099999 (1) w.o. free dual inf (0)
Clp0006I 0 Obj 0 Dual inf 0.0099999 (1) w.o. free dual inf (0)
Clp0006I 1 Obj -5
Clp0006I 1 Obj -5
Clp0000I Optimal - objective value -5
Clp0000I Optimal - objective value -5
Clp0006I 0 Obj 0 Dual inf 0.0149998 (2) w.o. free dual inf (0)
Clp0006I 0 Obj 0 Dual inf 0.0149998 (2) w.o. free dual inf (0)
Clp0006I 2 Obj -7.5
Clp0006I 2 Obj -7.5
Clp0000I Optimal - objective value -7.5
Clp0000I Optimal - objective value -7.5
Clp0006I 0 Obj 0 Dual inf 0.0184997 (3) w.o. free dual inf (0)
Clp0006I 0 Obj 0 Dual inf 0.0094997 (3) w.o. free dual inf (0)
Clp0006I 3 Obj -10
Clp0006I 3 Obj -5.5
Clp0000I Optimal - objective value -10
Clp0000I Optimal - objective value -5.5
Clp0006I 0 Obj 0 Dual inf 0.0149998 (2) w.o. free dual inf (0)
Clp0006I 0 Obj 0 Dual inf 0.0149998 (2) w.o. free dual inf (0)
Clp0006I 2 Obj -7.5
Clp0006I 2 Obj -7.5
Clp0000I Optimal - objective value -7.5
Clp0000I Optimal - objective value -7.5
Clp0006I 0 Obj 0 Dual inf 0.0184997 (3) w.o. free dual inf (0)
Clp0006I 0 Obj 0 Dual inf 0.0094997 (3) w.o. free dual inf (0)
Clp0006I 3 Obj -10
Clp0006I 3 Obj -5.5
Clp0000I Optimal - objective value -10
Clp0000I Optimal - objective value -5.5
We then solve the MPC problem and plot the results. The optimization problem will both optimize over the state and action trajectories, as well as the sequence of regions that the system will follow. For this reason, it is a mixed-integer optimization problem.
x_0 = [-3, 0]
sol_mixint = mpc.solve(pars={"x_0": x_0})
NLP0012I
Num Status Obj It time Location
NLP0014I 1 OPT 9 7 0.014244
NLP0014I 76 OPT 1150.7232 21 0.035136
NLP0012I
Num Status Obj It time Location
NLP0014I 1 OPT 1150.7232 13 0.021152
Cbc0012I Integer solution of 1150.7232 found by FPump after 0 iterations and 0 nodes (2.10 seconds)
NLP0014I 2 OPT 9 14 0.024257
NLP0014I 3 OPT 9 15 0.025437
NLP0014I 4 OPT 9 14 0.02406
NLP0014I 5 OPT 9 15 0.025301
NLP0014I 6 OPT 9 14 0.023863
NLP0014I 7 OPT 9 15 0.025659
NLP0014I 8 OPT 9 15 0.025372
NLP0014I 9 OPT 9 15 0.025566
NLP0014I 10 OPT 9 14 0.023906
NLP0014I 11 OPT 9 15 0.0257
NLP0014I 12 OPT 9 15 0.025567
NLP0014I 13 OPT 9 15 0.029773
NLP0014I 14 OPT 9 14 0.025729
NLP0014I 15 OPT 9 15 0.025645
NLP0014I 16 OPT 9 15 0.025791
NLP0014I 17 OPT 9 15 0.026166
NLP0014I 18 OPT 9 15 0.026053
NLP0014I 19 OPT 9 15 0.025361
NLP0014I 20 OPT 9 14 0.023944
NLP0012I
Num Status Obj It time Location
NLP0014I 21 OPT 9 15 0.025271
NLP0014I 22 OPT 9 14 0.024252
NLP0014I 23 OPT 9 15 0.025813
NLP0014I 24 OPT 9 15 0.025932
NLP0014I 25 OPT 9 15 0.026226
NLP0014I 26 OPT 9 14 0.024328
NLP0014I 27 OPT 9 15 0.025609
NLP0014I 28 OPT 9 15 0.026332
NLP0014I 29 OPT 9 15 0.026114
NLP0014I 30 OPT 9 15 0.026043
NLP0014I 31 OPT 9 15 0.025833
NLP0014I 32 OPT 9 15 0.025646
NLP0014I 33 OPT 9 15 0.026189
NLP0014I 34 OPT 9 16 0.027185
NLP0014I 35 OPT 9 15 0.028253
NLP0014I 36 OPT 9 16 0.028546
NLP0014I 37 OPT 9 14 0.024174
NLP0014I 38 INFEAS 0.6 11 0.022434
NLP0014I 39 OPT 26.977556 12 0.019651
NLP0014I 40 OPT 26.977556 12 0.020561
NLP0012I
Num Status Obj It time Location
NLP0014I 41 INFEAS 2.9994629 20 0.038062
NLP0014I 42 OPT 26.977556 10 0.019451
Cbc0010I After 0 nodes, 1 on tree, 1150.7232 best solution, best possible 26.977556 (3.16 seconds)
NLP0014I 43 OPT 26.977556 16 0.029339
NLP0014I 44 OPT 26.977556 16 0.029444
NLP0014I 45 INFEAS 0.59999999 21 0.040094
NLP0014I 46 OPT 92.437798 22 0.039507
NLP0014I 47 OPT 92.437796 24 0.044121
NLP0014I 48 OPT 92.437796 24 0.042967
NLP0014I 49 OPT 92.437796 23 0.041867
NLP0014I 50 OPT 92.437798 24 0.045443
NLP0014I 51 OPT 92.437796 22 0.041756
NLP0014I 52 OPT 92.437796 26 0.047082
NLP0014I 53 OPT 92.437796 23 0.041849
NLP0014I 54 OPT 92.437798 25 0.04678
NLP0014I 55 OPT 92.437798 2761 4.594424
NLP0014I 56 OPT 92.437797 26 0.04526
NLP0014I 57 OPT 92.437796 26 0.044605
NLP0014I 58 OPT 92.437798 28 0.050265
NLP0014I 59 INFEAS 0.16623375 43 0.079806
NLP0014I 60 INFEAS 0.17599999 43 0.077555
NLP0012I
Num Status Obj It time Location
NLP0014I 61 INFEAS 0.59999999 21 0.041997
NLP0014I 62 OPT 92.437796 28 0.051161
NLP0014I 63 OPT 92.437796 27 0.049118
NLP0014I 64 OPT 92.437796 24 0.045318
NLP0014I 65 OPT 92.437796 26 0.046727
NLP0014I 66 OPT 92.437796 26 0.047239
NLP0014I 67 OPT 92.437796 23 0.042809
NLP0014I 68 OPT 92.437798 25 0.044018
NLP0014I 69 OPT 92.437796 23 0.041232
NLP0014I 70 OPT 92.437796 24 0.043168
NLP0014I 71 OPT 92.437796 25 0.044317
NLP0014I 72 OPT 92.437796 24 0.042532
NLP0014I 73 OPT 92.437798 24 0.041627
NLP0014I 74 OPT 92.437797 28 0.048402
NLP0014I 75 INFEAS 0.19067846 51 0.091017
NLP0014I 76 INFEAS 0.25583871 44 0.077983
NLP0014I 77 OPT 92.437796 29 0.050563
NLP0014I 78 OPT 92.437798 27 0.047339
NLP0014I 79 OPT 92.437796 25 0.044027
NLP0014I 80 OPT 92.437797 26 0.045361
NLP0012I
Num Status Obj It time Location
NLP0014I 81 INFEAS 0.14814815 38 0.069936
NLP0014I 82 INFEAS 0.19067845 46 0.080724
NLP0014I 83 OPT 92.437796 26 0.045752
NLP0014I 84 OPT 92.437796 26 0.045788
NLP0014I 85 OPT 92.437798 25 0.043924
NLP0014I 86 OPT 92.437796 23 0.040298
NLP0014I 87 OPT 92.437798 26 0.04598
NLP0014I 88 OPT 92.437796 24 0.043582
NLP0014I 89 OPT 92.437798 27 0.046388
NLP0014I 90 OPT 92.437796 24 0.044184
NLP0014I 91 OPT 92.437796 25 0.044475
NLP0014I 92 OPT 92.437796 27 0.04674
NLP0014I 93 INFEAS 0.1228859 52 0.090966
NLP0014I 94 INFEAS 0.26666666 35 0.064554
NLP0014I 95 INFEAS 0.13647058 47 0.085764
NLP0014I 96 INFEAS 0.26666666 36 0.066167
NLP0014I 97 OPT 92.437796 24 0.041547
NLP0014I 98 OPT 92.437796 26 0.045763
NLP0014I 99 OPT 92.437796 24 0.042063
NLP0014I 100 OPT 92.437798 26 0.046012
NLP0012I
Num Status Obj It time Location
NLP0014I 101 OPT 92.437796 32 0.055183
NLP0014I 102 OPT 92.437796 24 0.041977
NLP0014I 103 OPT 92.437797 26 0.045741
NLP0014I 104 OPT 92.437796 28 0.048305
NLP0014I 105 INFEAS 0.2 40 0.072272
NLP0014I 106 INFEAS 0.11375707 46 0.081652
NLP0014I 107 OPT 92.437796 25 0.044513
NLP0014I 108 OPT 92.437796 24 0.043464
NLP0014I 109 OPT 92.437796 24 0.042463
NLP0014I 110 OPT 92.437796 25 0.044243
NLP0014I 111 OPT 92.437796 29 0.050506
NLP0014I 112 OPT 92.437798 44 0.075027
NLP0014I 113 OPT 92.437796 25 0.043163
NLP0014I 114 OPT 92.437796 24 0.042087
NLP0014I 115 INFEAS 0.16623375 38 0.068344
NLP0014I 116 INFEAS 0.17599999 42 0.075414
NLP0014I 117 INFEAS 0.16623375 37 0.067635
NLP0014I 118 INFEAS 0.17599998 48 0.086307
NLP0014I 119 OPT 92.437796 24 0.042745
NLP0014I 120 OPT 92.437796 25 0.043817
NLP0012I
Num Status Obj It time Location
NLP0014I 121 OPT 92.437796 25 0.044092
NLP0014I 122 OPT 92.437796 24 0.042466
NLP0014I 123 INFEAS 0.12851858 50 0.089465
NLP0014I 124 INFEAS 0.26666666 37 0.066847
NLP0014I 125 INFEAS 0.1285186 54 0.097541
NLP0014I 126 INFEAS 0.26666666 37 0.066463
NLP0014I 127 OPT 92.437797 27 0.047099
NLP0014I 128 OPT 92.437798 27 0.046198
NLP0014I 129 INFEAS 0.14814815 47 0.083859
NLP0014I 130 INFEAS 0.19067845 49 0.08896
NLP0014I 131 OPT 92.437798 26 0.046474
NLP0014I 132 OPT 92.437796 26 0.045223
NLP0014I 133 INFEAS 0.12851856 47 0.085095
NLP0014I 134 INFEAS 0.26666666 36 0.066037
NLP0014I 135 OPT 92.437796 26 0.045275
NLP0014I 136 OPT 92.437796 27 0.046308
NLP0014I 137 INFEAS 0.2 41 0.07467
NLP0014I 138 INFEAS 0.09509508 53 0.09468
NLP0014I 139 INFEAS 0.2 34 0.060643
NLP0014I 140 OPT 1150.7233 28 0.048881
NLP0012I
Num Status Obj It time Location
NLP0014I 141 OPT 92.437798 25 0.043536
NLP0014I 142 OPT 92.437796 24 0.041298
Cbc0010I After 100 nodes, 20 on tree, 1150.7232 best solution, best possible 92.437796 (13.15 seconds)
NLP0014I 143 INFEAS 0.23031914 50 0.097
NLP0014I 144 INFEAS 0.2268235 51 0.109778
NLP0014I 145 OPT 92.437796 26 0.048507
NLP0014I 146 OPT 92.437796 24 0.043776
NLP0014I 147 INFEAS 0.13333333 40 0.07345
NLP0014I 148 OPT 278.47937 27 0.047972
NLP0014I 149 OPT 1440.7231 30 0.053714
NLP0014I 150 INFEAS 0.17599998 40 0.073268
NLP0014I 151 OPT 92.437796 23 0.040464
NLP0014I 152 OPT 92.437796 25 0.043738
NLP0014I 153 OPT 92.437798 26 0.044739
NLP0014I 154 OPT 92.437796 25 0.043034
NLP0014I 155 INFEAS 0.13333333 40 0.071836
NLP0014I 156 OPT 1411.7664 37 0.063363
NLP0014I 157 OPT 92.437797 24 0.042441
NLP0014I 158 OPT 92.437798 27 0.047551
NLP0014I 159 INFEAS 0.16623375 36 0.064733
NLP0014I 160 INFEAS 0.17599998 42 0.075182
NLP0012I
Num Status Obj It time Location
NLP0014I 161 OPT 92.437796 25 0.043473
NLP0014I 162 OPT 92.437796 23 0.040573
NLP0014I 163 INFEAS 0.19613165 42 0.076061
NLP0014I 164 INFEAS 0.15797089 47 0.084486
NLP0014I 165 INFEAS 0.13333333 39 0.069893
NLP0014I 166 OPT 1062.5529 28 0.049456
NLP0014I 167 OPT 1064.9704 25 0.042476
NLP0014I 2 OPT 1064.9703 12 0.018494
Cbc0004I Integer solution of 1064.9703 found after 6582 iterations and 125 nodes (14.67 seconds)
NLP0014I 168 INFEAS 0.2 36 0.066143
NLP0014I 169 OPT 649.05551 32 0.055086
NLP0014I 170 INFEAS 0.13647058 45 0.081512
NLP0014I 171 OPT 1440.7231 39 0.068096
NLP0014I 172 INFEAS 0.2 43 0.077384
NLP0014I 173 INFEAS 0.11090709 50 0.089037
NLP0014I 174 INFEAS 0.14814815 47 0.084482
NLP0014I 175 INFEAS 0.19067846 51 0.091926
NLP0014I 176 INFEAS 0.14814815 44 0.078293
NLP0014I 177 INFEAS 0.19067845 51 0.089414
NLP0014I 178 INFEAS 0.16623375 41 0.075118
NLP0014I 179 INFEAS 0.17599999 43 0.076866
NLP0014I 180 OPT 92.437798 25 0.044142
NLP0012I
Num Status Obj It time Location
NLP0014I 181 OPT 92.437798 26 0.046782
NLP0014I 182 INFEAS 0.22710624 38 0.068717
NLP0014I 183 INFEAS 0.15210346 41 0.073405
NLP0014I 184 OPT 92.437796 32 0.055234
NLP0014I 185 OPT 92.437796 28 0.0473
NLP0014I 186 INFEAS 0.2 39 0.069398
NLP0014I 187 INFEAS 0.14631897 43 0.076813
NLP0014I 188 INFEAS 0.2 40 0.071613
NLP0014I 189 INFEAS 0.15189146 44 0.078746
NLP0014I 190 OPT 92.437796 23 0.040782
NLP0014I 191 OPT 92.437796 25 0.043968
NLP0014I 192 OPT 92.437796 23 0.041444
NLP0014I 193 OPT 92.437796 25 0.043494
NLP0014I 194 INFEAS 0.22985684 43 0.07789
NLP0014I 195 INFEAS 0.13647058 42 0.076381
NLP0014I 196 OPT 92.437796 24 0.043473
NLP0014I 197 OPT 92.437796 23 0.041414
NLP0014I 198 INFEAS 0.2 36 0.0663
NLP0014I 199 OPT 976.7376 28 0.049984
NLP0014I 200 OPT 1002.843 23 0.041973
NLP0012I
Num Status Obj It time Location
NLP0014I 201 OPT 1417.0833 29 0.050054
NLP0014I 202 INFEAS 0.11375707 46 0.082263
NLP0014I 203 INFEAS 0.15210346 42 0.088394
NLP0014I 204 INFEAS 0.2 35 0.065913
NLP0014I 205 OPT 727.86334 42 0.072258
NLP0014I 206 OPT 854.0901 27 0.047053
NLP0014I 207 OPT 798.77336 26 0.046449
NLP0014I 208 INFEAS 0.014620461 52 0.095185
NLP0014I 209 INFEAS 0.15210346 43 0.078779
NLP0014I 210 INFEAS 0.22985684 42 0.075992
NLP0014I 211 INFEAS 0.13647057 43 0.078954
NLP0014I 212 INFEAS 0.12851858 50 0.092274
NLP0014I 213 INFEAS 0.26666666 38 0.070389
NLP0014I 214 INFEAS 0.16623375 36 0.066412
NLP0014I 215 INFEAS 0.17599999 43 0.078
NLP0014I 216 OPT 92.437798 26 0.044857
NLP0014I 217 OPT 92.437798 26 0.043884
NLP0014I 218 INFEAS 0.16623375 39 0.07064
NLP0014I 219 INFEAS 0.17599999 42 0.074918
NLP0014I 220 OPT 92.437797 26 0.045023
NLP0012I
Num Status Obj It time Location
NLP0014I 221 OPT 92.437798 26 0.046315
NLP0014I 222 INFEAS 0.16623375 36 0.066079
NLP0014I 223 INFEAS 0.17599998 44 0.082799
NLP0014I 224 INFEAS 0.22985684 53 0.094655
NLP0014I 225 INFEAS 0.11428484 46 0.084063
NLP0014I 226 OPT 92.437798 26 0.047212
NLP0014I 227 OPT 92.437798 27 0.04828
NLP0014I 228 INFEAS 0.14814815 44 0.078463
NLP0014I 229 INFEAS 0.19067846 51 0.098821
NLP0014I 230 OPT 92.437797 27 0.047502
NLP0014I 231 OPT 92.437798 26 0.044613
NLP0014I 232 INFEAS 0.22985684 51 0.091978
NLP0014I 233 INFEAS 0.10881213 53 0.092687
NLP0014I 234 INFEAS 0.22985684 52 0.09193
NLP0014I 235 INFEAS 0.10881213 49 0.087037
NLP0014I 236 OPT 92.437796 29 0.049888
NLP0014I 237 OPT 92.437798 26 0.044033
NLP0014I 238 INFEAS 0.14814815 46 0.082727
NLP0014I 239 INFEAS 0.19067845 45 0.08226
NLP0014I 240 INFEAS 0.16623375 38 0.068713
NLP0012I
Num Status Obj It time Location
NLP0014I 241 INFEAS 0.17599999 42 0.075974
NLP0014I 242 INFEAS 0.14814815 43 0.078893
Cbc0010I After 200 nodes, 9 on tree, 1064.9703 best solution, best possible 92.437798 (19.79 seconds)
NLP0014I 243 INFEAS 0.19067845 49 0.089552
NLP0014I 244 INFEAS 0.16623375 39 0.073707
NLP0014I 245 INFEAS 0.17599998 40 0.074603
NLP0014I 246 INFEAS 0.14814815 44 0.079981
NLP0014I 247 INFEAS 0.19067846 53 0.094548
NLP0014I 248 INFEAS 0.19067845 51 0.091555
NLP0014I 249 INFEAS 0.25583871 43 0.079911
NLP0014I 250 INFEAS 0.22710624 37 0.067557
NLP0014I 251 INFEAS 0.15210346 41 0.07424
NLP0014I 252 OPT 92.437798 26 0.045255
NLP0014I 253 OPT 92.437796 24 0.042234
NLP0014I 254 INFEAS 0.16623375 38 0.068255
NLP0014I 255 INFEAS 0.17599998 45 0.082169
NLP0014I 256 INFEAS 0.16623375 41 0.076237
NLP0014I 257 INFEAS 0.17599998 44 0.079662
NLP0014I 258 INFEAS 0.16623375 42 0.076604
NLP0014I 259 INFEAS 0.17599998 46 0.082083
NLP0014I 260 OPT 92.437798 28 0.049255
NLP0012I
Num Status Obj It time Location
NLP0014I 261 OPT 92.437796 26 0.0456
NLP0014I 262 INFEAS 0.13834313 41 0.073634
NLP0014I 263 INFEAS 0.13419994 53 0.092571
NLP0014I 264 INFEAS 0.1691712 46 0.081476
NLP0014I 265 INFEAS 0.11760265 55 0.097998
NLP0014I 266 INFEAS 0.076752901 52 0.094225
NLP0014I 267 INFEAS 0.15210346 42 0.077591
NLP0014I 268 OPT 1560.5217 33 0.058188
Cbc0001I Search completed - best objective 1064.970288426228, took 10487 iterations and 226 nodes (21.76 seconds)
Cbc0032I Strong branching done 20 times (588 iterations), fathomed 0 nodes and fixed 2 variables
Cbc0035I Maximum depth 8, 0 variables fixed on reduced cost
CasADi - 2025-11-18 10:24:56 WARNING("solver_bonmin_Nlp12:nlp_grad failed: NaN detected for output grad_gamma_p, at (row 0, col 0).") [.../casadi/core/oracle_function.cpp:408]
CasADi - 2025-11-18 10:24:56 WARNING("Failed to calculate multipliers") [.../casadi/core/nlpsol.cpp:835]
solver_bonmin_Nlp12 : t_proc (avg) t_wall (avg) n_eval
nlp_f | 23.72ms ( 1.42us) 22.26ms ( 1.33us) 16731
nlp_g | 139.12ms ( 8.31us) 129.91ms ( 7.76us) 16732
nlp_grad_f | 22.71ms ( 2.03us) 20.42ms ( 1.82us) 11196
nlp_hess_l | 20.49ms ( 1.69us) 17.70ms ( 1.46us) 12149
nlp_jac_g | 119.34ms ( 9.21us) 118.02ms ( 9.11us) 12954
total | 21.96 s ( 21.96 s) 21.96 s ( 21.96 s) 1
Affine time-varying dynamics#
As stated above, when using csnlp.wrappers.PwaMpc.set_pwa_dynamics to specify
the PWA dynamics, the numerical solver will optimize also over the sequence of regions
that the system will follow, thus it must find the solution to a logical/integer
problem. This is often computationally expensive. But an alternative exists: to
specify a fixed switching sequence of regions manually/externally, and let the solver
only optimize the state-action trajectory. This is of course in general
computationally much cheaper.
csnlp.wrappers.PwaMpc also allows for defining the affine dynamics while
manually providing the sequence of regions the system should follow, rather than
letting the solver optimize it. The dynamics are thus time-varying affine. It is then
the user’s responsibility to specify reasonable switching sequences.
Building again the MPC, but this time, affine#
Now lets explore the setting in which the switching sequence is passed rather than
optimized. We build the MPC as before, but now using the
csnlp.wrappers.PwaMpc.set_affine_time_varying_dynamics method to set the
dynamics of the system instead. Note that, since the sequence is fixed, we do not need
a mixed-integer solver, but we can use any QP solver.
mpc = wrappers.PwaMpc(nlp=Nlp[cs.SX](sym_type="SX"), prediction_horizon=N)
x, _ = mpc.state("x", 2)
u, _ = mpc.action("u")
mpc.set_affine_time_varying_dynamics(pwa_system)
mpc.constraint("state_constraints", D1 @ x - E1, "<=", 0)
mpc.constraint("input_constraints", D2 @ u - E2, "<=", 0)
mpc.minimize(cs.sumsqr(x) + cs.sumsqr(u))
mpc.init_solver({"record_time": True}, "qrqp")
-------------------------------------------
This is casadi::QRQP
Number of variables: 32
Number of constraints: 96
Number of nonzeros in H: 32
Number of nonzeros in A: 176
Number of nonzeros in KKT: 480
Number of nonzeros in QR(V): 354
Number of nonzeros in QR(R): 758
We then set the switching sequence to be the optimal one (gathered from
the previous solution) via csnlp.wrappers.PwaMpc.set_switching_sequence, and
solve the ensuing QP problem for the same initial state.
Iter Sing fk |pr| con |du| var min_R con last_tau Note
0 0 0 3 32 5.1e-308 1 0.2 18 0
1 0 4.9e+02 3 96 2.6e-13 2 0.0022 48 0.97 Enforcing ubz, i=96
2 0 1e+03 0.2 92 5.7e-13 3 0.00058 48 1 Added ubz to reduce |pr|, i=92
3 1 1.1e+03 0.028 92 0.098 16 0.00058 48 0.86 Enforced ubz to reduce |du|, i=62
4 0 1.1e+03 0.028 92 0.098 16 0.00058 62 1 Dropped ubz for regularity, i=96
5 0 1.1e+03 0.028 92 0.098 16 0.0019 46 2.3e-05 Dropped ubz to reduce |du|, i=62
6 0 1.1e+03 4.8e-15 41 1.1e-13 4 0.0019 46 1 Converged
Effects of suboptimal sequences#
As aforementioned, the sequence now is specified by the user externally. This means that also suboptimal switching sequences can be passed. The solver will still find a solution, as long as the sequence is feasible, but the cost will be higher than when the optimal sequnce is passed or the sequence is part of the optimization.
subopt_sequence = opt_sequence.copy()
subopt_sequence[3] = 0
mpc.set_switching_sequence(subopt_sequence)
sol_qp_suboptimal = mpc.solve(pars={"x_0": x_0})
Iter Sing fk |pr| con |du| var min_R con last_tau Note
0 0 0 3 32 5.1e-308 1 0.17 18 0
1 0 4.1e+02 2.1 96 2.8e-13 2 0.0017 48 1 Added ubz to reduce |pr|, i=96
2 0 6.8e+02 1.1 57 7.4e-13 6 0.00072 57 1 Added ubz to reduce |pr|, i=57
3 0 1.2e+03 0.021 92 1.1e-12 6 0.00043 48 1 Added ubz to reduce |pr|, i=92
4 0 1.2e+03 1.1e-14 45 1.1e-13 0 0.00043 48 1 Converged
Results#
Let’s take a look at the optimality of the three solutions. Of course, we expect the mixed-integer solution to be the optimal one, the QP solution with the optimal sequence to be the same, and the QP solution with the suboptimal sequence to be worse.
print(f"Optimal mixed-integer cost: {sol_mixint.f}")
print(f"Optimal QP cost: {sol_qp.f}")
print(f"Suboptimal QP cost: {sol_qp_suboptimal.f}")
Optimal mixed-integer cost: 1064.9702884262279
Optimal QP cost: 1064.9704105314072
Suboptimal QP cost: 1150.7233547549013
However, we have gained some computational efficiency by not optimizing over the sequence of regions. This can be seen in the time taken to solve the problems.
print(f"Optimal mixed-integer time: {sol_mixint.stats['t_wall_total']}")
print(f"Optimal QP time: {sol_qp.stats['t_wall_total']}")
print(f"Suboptimal QP time: {sol_qp_suboptimal.stats['t_wall_total']}")
Optimal mixed-integer time: 21.96246314
Optimal QP time: 0.000376187
Suboptimal QP time: 0.000279145
We can also finally plot the three results (optimal mixed-integer, optimal QP, and suboptimal QP problem solutions).
_, axs = plt.subplots(1, 2, constrained_layout=True, sharey=True, figsize=(12, 5))
t = np.linspace(0, N, N + 1)
axs[0].step(t, sol_mixint.vals["x"].T, where="post")
axs[0].step(t[:-1], sol_mixint.vals["u"].T, where="post", color="C4")
axs[1].step(t, sol_qp.vals["x"].T, where="post")
axs[1].step(t, sol_qp_suboptimal.vals["x"].T, where="post", ls="--")
axs[1].step(t[:-1], sol_qp.vals["u"].T, where="post")
axs[1].step(t[:-1], sol_qp_suboptimal.vals["u"].T, where="post", ls="--")
axs[0].set_xlabel("Time step")
axs[0].set_title("Optimal mixed-integer solution")
axs[0].legend([r"$x^\text{MIQP}_1$", r"$x^\text{MIQP}_2$", r"$u^\text{MIQP}$"])
axs[0].set_xlabel("Time step")
axs[1].set_title("Optimal and suboptimal QP solutions")
axs[1].legend(
[
r"$x^\text{QP}_1$",
r"$x^\text{QP}_2$",
r"$u^\text{QP}$",
r"$x^\text{subQP}_1$",
r"$x^\text{subQP}_2$",
r"$u^\text{subQP}$",
]
)
plt.show()
Total running time of the script: (0 minutes 22.370 seconds)