.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "auto_examples/optimal_control/scenario_approach_for_mpc.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_scenario_approach_for_mpc.py: .. _scenario_based_mpc_example: Scenario-based MPC ================== In this example, we demonstrate how to build a scenario-based model predictive control (SCMPC) using :class:`csnlp.wrappers.ScenarioBasedMpc`. The problem is inspired by the numerical results of :cite:`schildbach_scenario_2014` shown in Tables 1 and 2, for :math:`R = R_1 = R_2 = 0`. .. GENERATED FROM PYTHON SOURCE LINES 14-15 We start with the usual imports. Let us also fix the RNG. .. GENERATED FROM PYTHON SOURCE LINES 15-27 .. code-block:: Python import casadi as cs import matplotlib.pyplot as plt import numpy as np import numpy.typing as npt from csnlp import Nlp from csnlp.wrappers import ScenarioBasedMpc np_random = np.random.default_rng(69) .. GENERATED FROM PYTHON SOURCE LINES 28-32 Setup ----- We then define the dynamics of the problem. The only peculiarity is that we have embedded a varying parameter in the dynamics, which is modelled as a disturbance. .. GENERATED FROM PYTHON SOURCE LINES 32-45 .. code-block:: Python def get_dynamics() -> cs.Function: x = cs.SX.sym("x", 2) u = cs.SX.sym("u", 2) d = cs.SX.sym("d", 3) # 2 disturbances + 1 varying parameter (modelled as dist.) theta = d[2] A = cs.blockcat([[0.7, -0.1 * (2 + theta)], [-0.1 * (3 + 2 * theta), 0.9]]) B = cs.DM.eye(2) x_next = A @ x + B @ u + d[:2] return cs.Function("F", [x, u, d], [x_next], ["x", "u", "d"], ["x+"]) .. GENERATED FROM PYTHON SOURCE LINES 46-48 In the Scenario Approach, we also need a function that allows us to sample independent samples of the disturbances affecting the stochastic optimization problem. .. GENERATED FROM PYTHON SOURCE LINES 48-56 .. code-block:: Python def sample_disturbances(K: int, N: int) -> npt.NDArray[np.floating]: d = np_random.normal(scale=0.1, size=(K, 2, N)) theta = np_random.uniform(size=(K, 1, N)) return np.concatenate((d, theta), 1) .. GENERATED FROM PYTHON SOURCE LINES 57-58 Let us also define some constants. .. GENERATED FROM PYTHON SOURCE LINES 58-64 .. code-block:: Python # parameters N = 5 # mpc horizon K = 19 # number of scenarios x0 = np.array([1, 1]) # initial state .. GENERATED FROM PYTHON SOURCE LINES 65-70 Building the SCMPC ------------------ The interface of the :class:`csnlp.wrappers.ScenarioBasedMpc` remains mostly similar to the one of :class:`csnlp.wrappers.Mpc`, so if this procedure is not clear, see the other, simpler example in :ref:`optimal_control_examples`. .. GENERATED FROM PYTHON SOURCE LINES 70-80 .. code-block:: Python F = get_dynamics() nx, nu, nd = F.size1_in(0), F.size1_in(1), F.size1_in(2) scmpc = ScenarioBasedMpc[cs.SX](Nlp(), K, N, shooting="single") x, _, _ = scmpc.state("x", nx, bound_initial=False) u, _ = scmpc.action("u", nu, lb=-5, ub=+5) d, _ = scmpc.disturbance("d", nd) scmpc.set_nonlinear_dynamics(F) .. GENERATED FROM PYTHON SOURCE LINES 81-87 Since the SCMPC internally defines ``K`` copies of the state and disturbances, how can we define a constraint that should be applied to each of the scenarios? We can use the method :meth:`csnlp.wrappers.ScenarioBasedMpc.constraint_from_single`, which will automatically translate the constraint to each scenario. In this case, we impose a lower bound on the state (we could have achieved the same result by passing `lb=1` to the state definition). .. GENERATED FROM PYTHON SOURCE LINES 87-90 .. code-block:: Python _ = scmpc.constraint_from_single("x_lb", x[:, 1:], ">=", 1) # equivalent to a lb on x .. GENERATED FROM PYTHON SOURCE LINES 91-95 The same reasoning is valid for the objective. It suffices to create the objective for one scenario, and then call :meth:`csnlp.wrappers.ScenarioBasedMpc.minimize_from_single`, which will automatically compute the average objective across all scenarios. .. GENERATED FROM PYTHON SOURCE LINES 95-100 .. code-block:: Python scmpc.minimize_from_single( sum(cs.sumsqr(x[:, i]) + cs.sumsqr(u[:, i]) for i in range(N)) ) .. GENERATED FROM PYTHON SOURCE LINES 101-104 Simulation ---------- We can finally simulate! First, let's initialize the solver. .. GENERATED FROM PYTHON SOURCE LINES 104-115 .. code-block:: Python opts = { "error_on_fail": True, "expand": True, "print_time": False, "record_time": True, "verbose": False, "printLevel": "none", } scmpc.init_solver(opts, "qpoases") .. rst-class:: sphx-glr-script-out .. code-block:: none qpOASES -- An Implementation of the Online Active Set Strategy. Copyright (C) 2007-2015 by Hans Joachim Ferreau, Andreas Potschka, Christian Kirches et al. All rights reserved. qpOASES is distributed under the terms of the GNU Lesser General Public License 2.1 in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more details. .. GENERATED FROM PYTHON SOURCE LINES 116-122 Then, we'd like to simulate the solution of the MPC along a trajectory of length ``T``. We'll save the resulting states and actions in dedicated lists. As commonly done, we adopt a receding horizon strategy, where only first action of the solution to the MPC is applied, and the rest is discarded. The process is then repeated at the next time step. .. GENERATED FROM PYTHON SOURCE LINES 122-145 .. code-block:: Python T = 200 # simulation repetitions x = x0 X, U = [], [] vals0 = None for t in range(T): # sample disturbances, and assign each to a scenario sample = sample_disturbances(scmpc.n_scenarios, scmpc.prediction_horizon) pars = {scmpc.name_i("d", i): sample[i] for i in range(scmpc.n_scenarios)} pars["x_0"] = x # set initial state # run the scenario-based MPC sol = scmpc.solve(pars, vals0) assert sol.success, f"Solver failed at time {t}!" u_opt = sol.vals["u"][:, 0] # apply only the first action # step the real system dynamics actual_disturbance_realization = sample_disturbances(1, 1).flatten() x = F(x, u_opt, actual_disturbance_realization).full().flatten() U.append(u_opt) X.append(x) .. rst-class:: sphx-glr-script-out .. code-block:: none qpOASES -- An Implementation of the Online Active Set Strategy. Copyright (C) 2007-2015 by Hans Joachim Ferreau, Andreas Potschka, Christian Kirches et al. All rights reserved. qpOASES is distributed under the terms of the GNU Lesser General Public License 2.1 in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more details. .. GENERATED FROM PYTHON SOURCE LINES 146-152 Results ------- First of all, we can compute the mean and standard deviation of the cost along the trajectory. The closer to zero, the better the performance of the controller. However, the controller should not be too aggressive, as it could lead to violating the constraints. .. GENERATED FROM PYTHON SOURCE LINES 152-157 .. code-block:: Python X, U = np.asarray(X), np.squeeze(U) Q, R = np.eye(nx), np.eye(nu) L = (X.dot(Q) * X).sum(1) + (U.dot(R) * U).sum(1) # same as the objective function print(f"Average cost: {L.mean():.2f} +/- {L.std():.2f}") .. rst-class:: sphx-glr-script-out .. code-block:: none Average cost: 3.83 +/- 0.57 .. GENERATED FROM PYTHON SOURCE LINES 158-162 We can then compute how many times either constraint (or both) has been violated. The higher the number of scenarios, the lower the chances of violations. In the plot, we show with marker ``x`` all the states that have violated a bound (different colors for different combinations of violations) .. GENERATED FROM PYTHON SOURCE LINES 162-179 .. code-block:: Python violated_con1, violated_con2 = X.T < 1 violated_either = np.logical_or(violated_con1, violated_con2) print(f"Average violation 1: {violated_con1.sum() / T * 100.0:.2f}%") print(f"Average violation 2: {violated_con2.sum() / T * 100.0:.2f}%") print(f"Average violation 1 or 2: {violated_either.sum() / T * 100.0:.2f}%") plt.axhline(1, color="darkgrey") plt.axvline(1, color="darkgrey") plt.plot(*X[~violated_con1 & ~violated_con2].T, "o", color="C0") plt.plot(*X[violated_con1 & ~violated_con2].T, "x", color="C1") plt.plot(*X[~violated_con1 & violated_con2].T, "x", color="C2") plt.plot(*X[violated_con1 & violated_con2].T, "x", color="C3") plt.axis("square") plt.xlabel("$x_1$") plt.ylabel("$x_2$") plt.show() .. image-sg:: /auto_examples/optimal_control/images/sphx_glr_scenario_approach_for_mpc_001.png :alt: scenario approach for mpc :srcset: /auto_examples/optimal_control/images/sphx_glr_scenario_approach_for_mpc_001.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none Average violation 1: 3.50% Average violation 2: 5.50% Average violation 1 or 2: 8.50% .. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 0.371 seconds) .. _sphx_glr_download_auto_examples_optimal_control_scenario_approach_for_mpc.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: scenario_approach_for_mpc.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: scenario_approach_for_mpc.py ` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: scenario_approach_for_mpc.zip ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_