.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "auto_examples/optimal_control/mpc_for_pwa_system.py" .. LINE NUMBERS ARE GIVEN BELOW. .. only:: html .. note:: :class: sphx-glr-download-link-note :ref:`Go to the end ` to download the full example code. .. rst-class:: sphx-glr-example-title .. _sphx_glr_auto_examples_optimal_control_mpc_for_pwa_system.py: MPC controller for PWA systems ============================== This example demos how an MPC controller can be built for a simple system with piecewise affine dynamics using :class:`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 :cite:`bemporad_control_1999` and :cite:`borrelli_predictive_2017`. 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 :math:`i`, the dynamics are described as .. math:: 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 .. math:: 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. .. GENERATED FROM PYTHON SOURCE LINES 30-32 We start with the imports as usual. All MPC controller classes can be found in the :mod:`csnlp.wrappers` module. .. GENERATED FROM PYTHON SOURCE LINES 32-39 .. code-block:: Python import casadi as cs import matplotlib.pyplot as plt import numpy as np from csnlp import Nlp, wrappers .. GENERATED FROM PYTHON SOURCE LINES 40-43 ----- Setup ----- .. GENERATED FROM PYTHON SOURCE LINES 45-50 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. .. GENERATED FROM PYTHON SOURCE LINES 50-60 .. code-block:: Python 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]]) .. GENERATED FROM PYTHON SOURCE LINES 61-63 From these matrices, we can define a sequence of :class:`csnlp.wrappers.PwaRegion` objects that represent the dynamics of the system in each region. .. GENERATED FROM PYTHON SOURCE LINES 63-85 .. code-block:: Python 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), ), ) .. GENERATED FROM PYTHON SOURCE LINES 86-91 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. .. GENERATED FROM PYTHON SOURCE LINES 91-101 .. code-block:: Python x_bnd = (5, 5) u_bnd = 20 D1 = np.array([[1, 0], [-1, 0], [0, 1], [0, -1]]) D2 = np.array([[1], [-1]]) D = cs.diagcat(D1, D2).sparse() E1 = np.array([x_bnd[0], x_bnd[0], x_bnd[1], x_bnd[1]]) E2 = np.array([u_bnd, u_bnd]) E = np.concatenate((E1, E2)) .. GENERATED FROM PYTHON SOURCE LINES 102-109 ------------------ 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 :meth:`csnlp.wrappers.PwaMpc.set_pwa_dynamics` instead of :meth:`csnlp.wrappers.Mpc.set_affine_dynamics` to set the dynamics of the system. .. GENERATED FROM PYTHON SOURCE LINES 109-120 .. code-block:: Python 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" .. rst-class:: sphx-glr-script-out .. code-block:: none 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 .. GENERATED FROM PYTHON SOURCE LINES 121-125 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. .. GENERATED FROM PYTHON SOURCE LINES 125-129 .. code-block:: Python x_0 = [-3, 0] sol_mixint = mpc.solve(pars={"x_0": x_0}) .. rst-class:: sphx-glr-script-out .. code-block:: none 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 .. GENERATED FROM PYTHON SOURCE LINES 130-144 ---------------------------- Affine time-varying dynamics ---------------------------- As stated above, when using :meth:`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. :class:`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. .. GENERATED FROM PYTHON SOURCE LINES 146-153 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 :meth:`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. .. GENERATED FROM PYTHON SOURCE LINES 153-163 .. code-block:: Python 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") .. rst-class:: sphx-glr-script-out .. code-block:: none ------------------------------------------- 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 .. GENERATED FROM PYTHON SOURCE LINES 164-167 We then set the switching sequence to be the optimal one (gathered from the previous solution) via :meth:`csnlp.wrappers.PwaMpc.set_switching_sequence`, and solve the ensuing QP problem for the same initial state. .. GENERATED FROM PYTHON SOURCE LINES 167-172 .. code-block:: Python opt_sequence = wrappers.PwaMpc.get_optimal_switching_sequence(sol_mixint) mpc.set_switching_sequence(opt_sequence) sol_qp = mpc.solve(pars={"x_0": x_0}) .. rst-class:: sphx-glr-script-out .. code-block:: none 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 .. GENERATED FROM PYTHON SOURCE LINES 173-179 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. .. GENERATED FROM PYTHON SOURCE LINES 179-185 .. code-block:: Python 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}) .. rst-class:: sphx-glr-script-out .. code-block:: none 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 .. GENERATED FROM PYTHON SOURCE LINES 186-192 ------- 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. .. GENERATED FROM PYTHON SOURCE LINES 192-197 .. code-block:: Python 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}") .. rst-class:: sphx-glr-script-out .. code-block:: none Optimal mixed-integer cost: 1064.9702884262279 Optimal QP cost: 1064.9704105314072 Suboptimal QP cost: 1150.7233547549013 .. GENERATED FROM PYTHON SOURCE LINES 198-200 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. .. GENERATED FROM PYTHON SOURCE LINES 200-205 .. code-block:: Python 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']}") .. rst-class:: sphx-glr-script-out .. code-block:: none Optimal mixed-integer time: 21.96246314 Optimal QP time: 0.000376187 Suboptimal QP time: 0.000279145 .. GENERATED FROM PYTHON SOURCE LINES 206-208 We can also finally plot the three results (optimal mixed-integer, optimal QP, and suboptimal QP problem solutions). .. GENERATED FROM PYTHON SOURCE LINES 208-236 .. code-block:: Python _, 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() .. image-sg:: /auto_examples/optimal_control/images/sphx_glr_mpc_for_pwa_system_001.png :alt: Optimal mixed-integer solution, Optimal and suboptimal QP solutions :srcset: /auto_examples/optimal_control/images/sphx_glr_mpc_for_pwa_system_001.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 22.370 seconds) .. _sphx_glr_download_auto_examples_optimal_control_mpc_for_pwa_system.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: mpc_for_pwa_system.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: mpc_for_pwa_system.py ` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: mpc_for_pwa_system.zip ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_