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.013913
NLP0014I 76 OPT 1150.7232 21 0.033251
NLP0012I
Num Status Obj It time Location
NLP0014I 1 OPT 1150.7232 13 0.02131
Cbc0012I Integer solution of 1150.7232 found by FPump after 0 iterations and 0 nodes (2.09 seconds)
NLP0014I 2 OPT 9 14 0.024895
NLP0014I 3 OPT 9 15 0.026256
NLP0014I 4 OPT 9 14 0.024877
NLP0014I 5 OPT 9 15 0.026299
NLP0014I 6 OPT 9 14 0.023985
NLP0014I 7 OPT 9 15 0.025743
NLP0014I 8 OPT 9 15 0.024584
NLP0014I 9 OPT 9 15 0.026046
NLP0014I 10 OPT 9 14 0.02415
NLP0014I 11 OPT 9 15 0.025676
NLP0014I 12 OPT 9 15 0.025789
NLP0014I 13 OPT 9 15 0.02586
NLP0014I 14 OPT 9 14 0.024182
NLP0014I 15 OPT 9 15 0.025719
NLP0014I 16 OPT 9 15 0.024649
NLP0014I 17 OPT 9 15 0.02582
NLP0014I 18 OPT 9 15 0.025912
NLP0014I 19 OPT 9 15 0.024743
NLP0014I 20 OPT 9 14 0.024221
NLP0012I
Num Status Obj It time Location
NLP0014I 21 OPT 9 15 0.025633
NLP0014I 22 OPT 9 14 0.02473
NLP0014I 23 OPT 9 15 0.02435
NLP0014I 24 OPT 9 15 0.026038
NLP0014I 25 OPT 9 15 0.025653
NLP0014I 26 OPT 9 14 0.024886
NLP0014I 27 OPT 9 15 0.025991
NLP0014I 28 OPT 9 15 0.025603
NLP0014I 29 OPT 9 15 0.025748
NLP0014I 30 OPT 9 15 0.025773
NLP0014I 31 OPT 9 15 0.025611
NLP0014I 32 OPT 9 15 0.026225
NLP0014I 33 OPT 9 15 0.025235
NLP0014I 34 OPT 9 16 0.027117
NLP0014I 35 OPT 9 15 0.027998
NLP0014I 36 OPT 9 16 0.027444
NLP0014I 37 OPT 9 14 0.023926
NLP0014I 38 INFEAS 0.6 11 0.021598
NLP0014I 39 OPT 26.977556 12 0.020395
NLP0014I 40 OPT 26.977556 12 0.019493
NLP0012I
Num Status Obj It time Location
NLP0014I 41 INFEAS 2.9994629 20 0.037493
NLP0014I 42 OPT 26.977556 10 0.019739
Cbc0010I After 0 nodes, 1 on tree, 1150.7232 best solution, best possible 26.977556 (3.14 seconds)
NLP0014I 43 OPT 26.977556 16 0.029967
NLP0014I 44 OPT 26.977556 16 0.029563
NLP0014I 45 INFEAS 0.59999999 21 0.041016
NLP0014I 46 OPT 92.437798 22 0.039321
NLP0014I 47 OPT 92.437796 24 0.042956
NLP0014I 48 OPT 92.437796 24 0.042947
NLP0014I 49 OPT 92.437796 23 0.042262
NLP0014I 50 OPT 92.437798 24 0.043522
NLP0014I 51 OPT 92.437796 22 0.039855
NLP0014I 52 OPT 92.437796 26 0.046028
NLP0014I 53 OPT 92.437796 23 0.040423
NLP0014I 54 OPT 92.437798 25 0.046521
NLP0014I 55 OPT 92.437798 2761 4.607461
NLP0014I 56 OPT 92.437797 26 0.049764
NLP0014I 57 OPT 92.437796 26 0.045095
NLP0014I 58 OPT 92.437798 28 0.048323
NLP0014I 59 INFEAS 0.16623375 43 0.07752
NLP0014I 60 INFEAS 0.17599999 43 0.077208
NLP0012I
Num Status Obj It time Location
NLP0014I 61 INFEAS 0.59999999 21 0.041834
NLP0014I 62 OPT 92.437796 28 0.050719
NLP0014I 63 OPT 92.437796 27 0.047907
NLP0014I 64 OPT 92.437796 24 0.044208
NLP0014I 65 OPT 92.437796 26 0.045191
NLP0014I 66 OPT 92.437796 26 0.044927
NLP0014I 67 OPT 92.437796 23 0.041077
NLP0014I 68 OPT 92.437798 25 0.043573
NLP0014I 69 OPT 92.437796 23 0.041216
NLP0014I 70 OPT 92.437796 24 0.042666
NLP0014I 71 OPT 92.437796 25 0.044164
NLP0014I 72 OPT 92.437796 24 0.042249
NLP0014I 73 OPT 92.437798 24 0.04052
NLP0014I 74 OPT 92.437797 28 0.048403
NLP0014I 75 INFEAS 0.19067846 51 0.090403
NLP0014I 76 INFEAS 0.25583871 44 0.077791
NLP0014I 77 OPT 92.437796 29 0.050747
NLP0014I 78 OPT 92.437798 27 0.045787
NLP0014I 79 OPT 92.437796 25 0.042331
NLP0014I 80 OPT 92.437797 26 0.044889
NLP0012I
Num Status Obj It time Location
NLP0014I 81 INFEAS 0.14814815 38 0.068546
NLP0014I 82 INFEAS 0.19067845 46 0.082325
NLP0014I 83 OPT 92.437796 26 0.046489
NLP0014I 84 OPT 92.437796 26 0.046231
NLP0014I 85 OPT 92.437798 25 0.044271
NLP0014I 86 OPT 92.437796 23 0.040135
NLP0014I 87 OPT 92.437798 26 0.045346
NLP0014I 88 OPT 92.437796 24 0.042427
NLP0014I 89 OPT 92.437798 27 0.047608
NLP0014I 90 OPT 92.437796 24 0.043092
NLP0014I 91 OPT 92.437796 25 0.042712
NLP0014I 92 OPT 92.437796 27 0.048105
NLP0014I 93 INFEAS 0.1228859 52 0.092123
NLP0014I 94 INFEAS 0.26666666 35 0.064119
NLP0014I 95 INFEAS 0.13647058 47 0.085313
NLP0014I 96 INFEAS 0.26666666 36 0.064597
NLP0014I 97 OPT 92.437796 24 0.042118
NLP0014I 98 OPT 92.437796 26 0.04583
NLP0014I 99 OPT 92.437796 24 0.042976
NLP0014I 100 OPT 92.437798 26 0.046989
NLP0012I
Num Status Obj It time Location
NLP0014I 101 OPT 92.437796 32 0.056589
NLP0014I 102 OPT 92.437796 24 0.042655
NLP0014I 103 OPT 92.437797 26 0.045654
NLP0014I 104 OPT 92.437796 28 0.048072
NLP0014I 105 INFEAS 0.2 40 0.071692
NLP0014I 106 INFEAS 0.11375707 46 0.082899
NLP0014I 107 OPT 92.437796 25 0.044774
NLP0014I 108 OPT 92.437796 24 0.041778
NLP0014I 109 OPT 92.437796 24 0.042771
NLP0014I 110 OPT 92.437796 25 0.044515
NLP0014I 111 OPT 92.437796 29 0.050445
NLP0014I 112 OPT 92.437798 44 0.073813
NLP0014I 113 OPT 92.437796 25 0.043697
NLP0014I 114 OPT 92.437796 24 0.042378
NLP0014I 115 INFEAS 0.16623375 38 0.070179
NLP0014I 116 INFEAS 0.17599999 42 0.075285
NLP0014I 117 INFEAS 0.16623375 37 0.06779
NLP0014I 118 INFEAS 0.17599998 48 0.087305
NLP0014I 119 OPT 92.437796 24 0.043909
NLP0014I 120 OPT 92.437796 25 0.045133
NLP0012I
Num Status Obj It time Location
NLP0014I 121 OPT 92.437796 25 0.04494
NLP0014I 122 OPT 92.437796 24 0.042181
NLP0014I 123 INFEAS 0.12851858 50 0.090954
NLP0014I 124 INFEAS 0.26666666 37 0.067296
NLP0014I 125 INFEAS 0.1285186 54 0.09784
NLP0014I 126 INFEAS 0.26666666 37 0.067745
NLP0014I 127 OPT 92.437797 27 0.046267
NLP0014I 128 OPT 92.437798 27 0.04676
NLP0014I 129 INFEAS 0.14814815 47 0.084689
NLP0014I 130 INFEAS 0.19067845 49 0.087978
NLP0014I 131 OPT 92.437798 26 0.044476
NLP0014I 132 OPT 92.437796 26 0.046493
NLP0014I 133 INFEAS 0.12851856 47 0.086117
NLP0014I 134 INFEAS 0.26666666 36 0.06757
NLP0014I 135 OPT 92.437796 26 0.045212
NLP0014I 136 OPT 92.437796 27 0.047021
NLP0014I 137 INFEAS 0.2 41 0.074508
NLP0014I 138 INFEAS 0.09509508 53 0.094459
NLP0014I 139 INFEAS 0.2 34 0.061758
NLP0014I 140 OPT 1150.7233 28 0.049041
NLP0012I
Num Status Obj It time Location
NLP0014I 141 OPT 92.437798 25 0.043542
NLP0014I 142 OPT 92.437796 24 0.04268
Cbc0010I After 100 nodes, 20 on tree, 1150.7232 best solution, best possible 92.437796 (13.13 seconds)
NLP0014I 143 INFEAS 0.23031914 50 0.089658
NLP0014I 144 INFEAS 0.2268235 51 0.091411
NLP0014I 145 OPT 92.437796 26 0.05258
NLP0014I 146 OPT 92.437796 24 0.05073
NLP0014I 147 INFEAS 0.13333333 40 0.088275
NLP0014I 148 OPT 278.47937 27 0.047349
NLP0014I 149 OPT 1440.7231 30 0.052352
NLP0014I 150 INFEAS 0.17599998 40 0.073309
NLP0014I 151 OPT 92.437796 23 0.041463
NLP0014I 152 OPT 92.437796 25 0.044459
NLP0014I 153 OPT 92.437798 26 0.044032
NLP0014I 154 OPT 92.437796 25 0.043289
NLP0014I 155 INFEAS 0.13333333 40 0.071898
NLP0014I 156 OPT 1411.7664 37 0.062372
NLP0014I 157 OPT 92.437797 24 0.042359
NLP0014I 158 OPT 92.437798 27 0.048192
NLP0014I 159 INFEAS 0.16623375 36 0.066193
NLP0014I 160 INFEAS 0.17599998 42 0.076111
NLP0012I
Num Status Obj It time Location
NLP0014I 161 OPT 92.437796 25 0.043898
NLP0014I 162 OPT 92.437796 23 0.041258
NLP0014I 163 INFEAS 0.19613165 42 0.076064
NLP0014I 164 INFEAS 0.15797089 47 0.083379
NLP0014I 165 INFEAS 0.13333333 39 0.072586
NLP0014I 166 OPT 1062.5529 28 0.050593
NLP0014I 167 OPT 1064.9704 25 0.043937
NLP0014I 2 OPT 1064.9703 12 0.020959
Cbc0004I Integer solution of 1064.9703 found after 6582 iterations and 125 nodes (14.67 seconds)
NLP0014I 168 INFEAS 0.2 36 0.066362
NLP0014I 169 OPT 649.05551 32 0.055021
NLP0014I 170 INFEAS 0.13647058 45 0.083259
NLP0014I 171 OPT 1440.7231 39 0.070267
NLP0014I 172 INFEAS 0.2 43 0.077743
NLP0014I 173 INFEAS 0.11090709 50 0.089249
NLP0014I 174 INFEAS 0.14814815 47 0.083671
NLP0014I 175 INFEAS 0.19067846 51 0.0921
NLP0014I 176 INFEAS 0.14814815 44 0.078815
NLP0014I 177 INFEAS 0.19067845 51 0.090471
NLP0014I 178 INFEAS 0.16623375 41 0.074892
NLP0014I 179 INFEAS 0.17599999 43 0.07612
NLP0014I 180 OPT 92.437798 25 0.044414
NLP0012I
Num Status Obj It time Location
NLP0014I 181 OPT 92.437798 26 0.048142
NLP0014I 182 INFEAS 0.22710624 38 0.071591
NLP0014I 183 INFEAS 0.15210346 41 0.076046
NLP0014I 184 OPT 92.437796 32 0.056062
NLP0014I 185 OPT 92.437796 28 0.048508
NLP0014I 186 INFEAS 0.2 39 0.06986
NLP0014I 187 INFEAS 0.14631897 43 0.07881
NLP0014I 188 INFEAS 0.2 40 0.071059
NLP0014I 189 INFEAS 0.15189146 44 0.078437
NLP0014I 190 OPT 92.437796 23 0.039695
NLP0014I 191 OPT 92.437796 25 0.044822
NLP0014I 192 OPT 92.437796 23 0.040928
NLP0014I 193 OPT 92.437796 25 0.043886
NLP0014I 194 INFEAS 0.22985684 43 0.077401
NLP0014I 195 INFEAS 0.13647058 42 0.078102
NLP0014I 196 OPT 92.437796 24 0.043341
NLP0014I 197 OPT 92.437796 23 0.040993
NLP0014I 198 INFEAS 0.2 36 0.063875
NLP0014I 199 OPT 976.7376 28 0.049151
NLP0014I 200 OPT 1002.843 23 0.041846
NLP0012I
Num Status Obj It time Location
NLP0014I 201 OPT 1417.0833 29 0.050053
NLP0014I 202 INFEAS 0.11375707 46 0.082481
NLP0014I 203 INFEAS 0.15210346 42 0.076411
NLP0014I 204 INFEAS 0.2 35 0.064488
NLP0014I 205 OPT 727.86334 42 0.072256
NLP0014I 206 OPT 854.0901 27 0.04826
NLP0014I 207 OPT 798.77336 26 0.047209
NLP0014I 208 INFEAS 0.014620461 52 0.094098
NLP0014I 209 INFEAS 0.15210346 43 0.080039
NLP0014I 210 INFEAS 0.22985684 42 0.076786
NLP0014I 211 INFEAS 0.13647057 43 0.07824
NLP0014I 212 INFEAS 0.12851858 50 0.089757
NLP0014I 213 INFEAS 0.26666666 38 0.069734
NLP0014I 214 INFEAS 0.16623375 36 0.066897
NLP0014I 215 INFEAS 0.17599999 43 0.077921
NLP0014I 216 OPT 92.437798 26 0.044979
NLP0014I 217 OPT 92.437798 26 0.045996
NLP0014I 218 INFEAS 0.16623375 39 0.07068
NLP0014I 219 INFEAS 0.17599999 42 0.07522
NLP0014I 220 OPT 92.437797 26 0.045237
NLP0012I
Num Status Obj It time Location
NLP0014I 221 OPT 92.437798 26 0.045375
NLP0014I 222 INFEAS 0.16623375 36 0.066682
NLP0014I 223 INFEAS 0.17599998 44 0.078752
NLP0014I 224 INFEAS 0.22985684 53 0.09482
NLP0014I 225 INFEAS 0.11428484 46 0.084845
NLP0014I 226 OPT 92.437798 26 0.046053
NLP0014I 227 OPT 92.437798 27 0.046553
NLP0014I 228 INFEAS 0.14814815 44 0.078125
NLP0014I 229 INFEAS 0.19067846 51 0.090749
NLP0014I 230 OPT 92.437797 27 0.047981
NLP0014I 231 OPT 92.437798 26 0.048244
NLP0014I 232 INFEAS 0.22985684 51 0.091382
NLP0014I 233 INFEAS 0.10881213 53 0.093736
NLP0014I 234 INFEAS 0.22985684 52 0.091103
NLP0014I 235 INFEAS 0.10881213 49 0.092589
NLP0014I 236 OPT 92.437796 29 0.049317
NLP0014I 237 OPT 92.437798 26 0.043861
NLP0014I 238 INFEAS 0.14814815 46 0.081303
NLP0014I 239 INFEAS 0.19067845 45 0.080614
NLP0014I 240 INFEAS 0.16623375 38 0.069005
NLP0012I
Num Status Obj It time Location
NLP0014I 241 INFEAS 0.17599999 42 0.075489
NLP0014I 242 INFEAS 0.14814815 43 0.078737
Cbc0010I After 200 nodes, 9 on tree, 1064.9703 best solution, best possible 92.437798 (19.78 seconds)
NLP0014I 243 INFEAS 0.19067845 49 0.088144
NLP0014I 244 INFEAS 0.16623375 39 0.070061
NLP0014I 245 INFEAS 0.17599998 40 0.073089
NLP0014I 246 INFEAS 0.14814815 44 0.079987
NLP0014I 247 INFEAS 0.19067846 53 0.093985
NLP0014I 248 INFEAS 0.19067845 51 0.092001
NLP0014I 249 INFEAS 0.25583871 43 0.07653
NLP0014I 250 INFEAS 0.22710624 37 0.069617
NLP0014I 251 INFEAS 0.15210346 41 0.077885
NLP0014I 252 OPT 92.437798 26 0.046426
NLP0014I 253 OPT 92.437796 24 0.041164
NLP0014I 254 INFEAS 0.16623375 38 0.07004
NLP0014I 255 INFEAS 0.17599998 45 0.081185
NLP0014I 256 INFEAS 0.16623375 41 0.07317
NLP0014I 257 INFEAS 0.17599998 44 0.080188
NLP0014I 258 INFEAS 0.16623375 42 0.074373
NLP0014I 259 INFEAS 0.17599998 46 0.080998
NLP0014I 260 OPT 92.437798 28 0.048272
NLP0012I
Num Status Obj It time Location
NLP0014I 261 OPT 92.437796 26 0.045113
NLP0014I 262 INFEAS 0.13834313 41 0.07329
NLP0014I 263 INFEAS 0.13419994 53 0.094186
NLP0014I 264 INFEAS 0.1691712 46 0.08103
NLP0014I 265 INFEAS 0.11760265 55 0.095111
NLP0014I 266 INFEAS 0.076752901 52 0.092712
NLP0014I 267 INFEAS 0.15210346 42 0.078307
NLP0014I 268 OPT 1560.5217 33 0.057818
Cbc0001I Search completed - best objective 1064.970288426228, took 10487 iterations and 226 nodes (21.73 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-10-16 20:31:52 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-10-16 20:31:52 WARNING("Failed to calculate multipliers") [.../casadi/core/nlpsol.cpp:835]
solver_bonmin_Nlp12 : t_proc (avg) t_wall (avg) n_eval
nlp_f | 25.55ms ( 1.53us) 24.01ms ( 1.43us) 16731
nlp_g | 142.90ms ( 8.54us) 134.24ms ( 8.02us) 16732
nlp_grad_f | 22.27ms ( 1.99us) 20.28ms ( 1.81us) 11196
nlp_hess_l | 21.37ms ( 1.76us) 19.32ms ( 1.59us) 12149
nlp_jac_g | 122.79ms ( 9.48us) 121.55ms ( 9.38us) 12954
total | 22.00 s ( 22.00 s) 22.00 s ( 22.00 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.999647317
Optimal QP time: 0.000374096
Suboptimal QP time: 0.000331376
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.410 seconds)