.. 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.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 .. 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.999647317 Optimal QP time: 0.000374096 Suboptimal QP time: 0.000331376 .. 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.410 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 `_