Welcome to nMPyC’s documentation!

nMPyC is a Python library for solving optimal control problems via model predictive control (MPC).

nMPyC can be understood as a blackbox method. The user can only enter the desired optimal control problem without having much knowledge of the theory of model predictive control or its implementation in Python. Nevertheless, for an advanced user, there is the possibility to adjust all parameters.

This library supports a variety of discretization methods and optimizers including CasADi and SciPy solvers.

In summary, nMPyC

solves nonlinear finite horizon optimal control problems
solves nonlinear optimal control problems with model predicitve control (MPC)
uses algorithmic differentation via CasADi
can chose between different discretization methods
can chose between different solvers for nonlinear optimization (depending on the problem)
supports time-varying optimal control problems
supports the special structure of linear-quadratic optimal control problems
supports discounted optimal control problems
can save and load the simulation results

The nMPyC software is Python based and works therefore on any OS with a Python distribution (for more precise requiremnents see the Installation section). nMPyC has been developed by Jonas Schießl and Lisa Krügel under the supervision of Prof. Lars Grüne at the Chair of Applied Mathematic of University of Bayreuth. nMPyC is a further devolpement in Python of the Matlab code that was implemented for the NMPC Book from Lars Grüne and Jürgen Pannek [GruneP17].

Contents

Installation

Requirements

The nMPyC package is dependent on the following libraries:

Installation using PIP

The easiest way to install the nMPyC package is to use PIP. To do so, you just need to run the command line

pip install nmpyc

The main advantage of this method is that the package is automatically added to the Python default path and all dependencies are installed.

Additionally you can update the package by running

pip install nmpyc --upgrade

Installation by Source

To install the Python package by source, the source code from GitHub has to be downloaded. This can be done via Git using the command

git clone https://github.com/nMPyC/nmpyc

Now the toolbox can be used by importing the package according to its storage path in the Python code by adding it to the Python default path. To realize the letter case you can navigate to the location of the package and use

pip install .

This command will automatically add the package to the Python default path and install the required Python packages and their dependencies.

Getting Started

Import nMPyC

After the successfull installation of the nMPyC package, nMPyC has to be imported to our code. This can be done as shown in the following code snippet

# Add nMPyC to path (if necessary)
import sys
sys.path.append('../../path-to-nmpyc')

# Import nmpyc
import nmpyc

Note that the first two lines can be omitted if nMPyC has already been added to the Python default path as described in the Installation section. In this case the command import nmpyc is sufficient to import the nMPyC library.

Note

Please use the nmpyc.nmpyc_array functions and the nmpyc.nmpyc_array.array class for the calculations in the code to ensure error-free functionality of the program. Further informations about this issue can be found in API References and in the FAQ section.

Creating the System Dynamics

To define the system dynamics of the optimal control problem, we have to create a nmpyc.system object. We can define the possibly time-dependent and nonlinear system dynamics using a function of the form

def f(t,x,u):
    y = nmpyc.array(nx)
    ...
    return y

If this function is created, the system can be initialized by calling

system = nmpyc.system(f,nx,nu,system_type)

Where nx is the dimension of the state, nu is the dimension of the control variable, and system_type is a string indicating whether the system is continuous (continuous) or discrete (discrete).

Furthermore, the parameters sampling_rate (sampling rate), t0 (initial time) and method can optionally be adjusted during the initialization of the system. The value of method determines the used integration method for the discretization of the differential equation in the continuous case. By default the CasADi integrator cvodes is used.

Further options of the used integration method can be defined by the command

system.set_integratorOptions(dict())

For more informations (also about the parameters and their standard values) see the API-References nmpyc.system.system.

Creating the Objective

To define the objective, we need to create – similar to the system dynamics – a nmpyc.objective object. To do so, we first define the stage cost

def l(t,x,u):
    ...
    return y

and add, optionally, a terminal cost of the form

def F(t,x):
    ...
    return y

Now we can initialize the objective by calling

objective = nmpyc.objective(l, F)
# Or alternatively without terminal costs
objective = nmpyc.objective(l)

For more informations see the API-References nmpyc.objective.objective.

Creating the Constraints

The optimal control problem can be extended with other constraints besides the necessary system dynamics. For this reason, we must first create an empty nmpyc.constraints object using the command

system = nmpyc.constraints()

We can now add the desired constraints to this object step by step. These constraints can be created in different ways. First, we can add box constraints in the form of bounds.

constraints.add_bound('lower', 'control', lbu) # lower bound for control
constraints.add_bound('upper', 'control', ubu) # upper bound for control

Here lbu or lbx is an nmpyc.nmpyc_array.array of dimension (1,nu) or (nu,1). To add bounds for the state or terminal state, replace control with state or terminal in the above code and adjust the dimension of the array accordingly.

In addition to box constraints, general inequality and equality constraints can also be inserted.

# Equality constraint h(t,x,u) = 0
def h(t,x,u):
   y = mpc.array(len_constr)
   ...
   return y
constraints.add_constr('eq', h)

# Inequality constraint g(t,x,u) >= 0
def g(t,x,u):
   y = mpc.array(len_constr)
   ...
   return y
constraints.add_constr('ineq', g)

Terminal constraints of the form \(H(t,x) = 0\) or \(G(t,x) \geq 0\) can also be added.

constraints.add_constr('terminal_eq', H)
constraints.add_constr('terminal_ineq', G)

Moreover it is possible to add linear equality and inequality constraints. For this purpose see nmpyc.constraints.constraints.add_constr(). For further general informations see the API-References nmpyc.constraints.constraints.

Running Simulations

After initializing all necessary objects, we can run simulations for our problem. We first create a mpc.model object and combine the different parts of the optimal control problem by calling

model = nmpyc.model(objective, system, constraints)

The nmpyc.constraints object is optional and can be omitted for a problem without constraints. Modyfying the default settings of the optimization, can be done with the help of the commands

model.opti.set_options(dict())
model.opti.set_solverOptions(dict())

For more informations about this methods see nmpyc.model.model.opti.

To start an open loop simulation, we execute the command

u_ol, x_ol = model.solve_ocp(x0,N,discount)

and for a closed loop simulation

res = model.mpc(x0,N,K,discount)

Here x0 is a nmpyc.nmpyc_array.array which defines the initial value, N is the MPC horizon and the parameter K defines the number of MPC iterations. The parameter discount is optional and defines the discount factor (the default is 1).

The result of the simulation can now be shown in the console by calling

print(res)

and as a visual output by calling

res.plot()

By default, the states and controls are displayed in two subplots. By passing a string as the first parameter (=args), the plot can be customized. For example, by calling

res.plot('state')

only the states are plotted. Other keywords are control for the control, cost for the stage costs, and phase to make a phase portrait of two states or controls. Furthermore, the plots displayed in this way can be additionally adjusted by further prameters, see nmpyc.result.result.plot().

Further, the model and the simulation results can be saved for later use with the functions

model.save('path')
res.save('path')

These saved files can then be loaded with the help of

model = nmpyc.model.load('path')
res = nmpyc.result.load('path')

Advanced topics

The procedure described above is only an excerpt of the possibilities of the nMPyC Python library. For example, it is also possible to create autonomous systems and use the linear quadratic structure of a problem. For further informations see the Examples and Templates section. And for the implementation of linear system dynamics and quadratic costs, see also nmpyc.system.system.LQP() and nmpyc.objective.objective.LQP().

Basics of model predictive control

Model predictive control (MPC) is an optimized-based method for obtaining an approximately optimal feedback control for an optimal control problem on an infinite or finite time horizon. The basic idea of MPC is to predict the future behavior of the controlled system over a finite time horizon and compute an optimal control input that, while ensuring satisfaction of given system constraints, minimizes the objective function. In each sampling instant a finite horizon open-loop optimal control problem is solved to calculate the control input. More precisley, this control input is used to define the feedback which is applied to the system until the next sampling instant, at whicht the horizon is shifted and the procedure is repeated again.

Optimal control problems

In order to describe the functionality of MPC we consider optimal control problems. To this end, we consider possibly nonlinear difference equations of the form

\[\begin{split}x(k+1,x_0) &= f(x(k,x_0),u(k)), \quad k = 0,\dots,N-1, \\ x(0) &= x_0\end{split}\]

with \(N\in\mathbb{N}\) or discretized differential equations.

Further, we impose nonempty state and input constraint sets \(\mathbb{X}\subseteq\mathbb{R}^{n}\) and \(\mathbb{U}\subseteq\mathbb{R}^m\), respectively, as well as a nonempty terminal constraint set \(\mathbb{X}_0\subseteq\mathbb{R}^n\).

Now we use optimal control to determine \(u(0),\dots,u(N-1)\). For this reason, we fix a stage cost \(\ell:\mathbb{X}\times\mathbb{U}\to\mathbb{R}\) which may be a very general function and a optional terminal cost \(F:\mathbb{X}\to\mathbb{R}\). Regardless which cost function is used the objective function is defined by

\[J^N(x_0,u(\cdot)):=\sum_{k=0}^{N-1}\ell(x(k,x_0),u(k))\]

without terminal cost or by

\[J^N(x_0,u(\cdot)):=\sum_{k=0}^{N-1}\ell(x(k,x_0),u(k))+ F(x(N,x_0))\]

with terminal cost.

In summary, an optimal control problem without terminal conditions is given by

\begin{equation} \begin{split} \min_{u(\cdot)\in\mathbb{U}}J^N(x_0,u(\cdot)) &= \sum_{k=0}^{N-1}\ell(x(k,x_0),u(k))\\ \text{s.t.}\quad x(k+1,x_0)&=f(x(k,x_0),u(k)),\quad k = 0,\dots, N-1\\ x(0)&= x_0\\ x&\in\mathbb{X} \end{split} \end{equation}

and an optimal control problem with terminal conditions is given by

\begin{equation} \begin{split} \min_{u(\cdot)\in\mathbb{U}}J^N(x_0,u(\cdot)) &= \sum_{k=0}^{N-1}\ell(x(k,x_0),u(k))+F(x(N,x_0))\\ \text{s.t.}\quad x(k+1,x_0)&=f(x(k,x_0),u(k)),\quad k = 0,\dots, N-1\\ x(0)&= x_0\\ x\in\mathbb{X},\quad & x(N,x_0)\in\mathbb{X}_0 \end{split} \end{equation}

Additionally, with nMPyC it is possible to add constraints to the optimal control problem.

The basic MPC algorithm

Regardless of the type of the optimal control problem, the MPC algorithm is given by:

At each time instant \(j=0,1,2,\dots:\)

Measure the state \(x(j)\in\mathbb{X}\) of the system.
Set \(x_0:=x(j)\), solve the optimal control problem (with or without terminal conditions) and denote the obtained optimal control sequence by \(u^\star(\cdot)\in\mathbb{U}^N(x_0)\).
Define the MPC-feedback value \(\mu^N(x(j)):=u^\star(0)\in\mathbb{U}\) and use this control value in the next sampling period (apply the feedback to the system).

Notes and extensions

A special case of an optimal control problem is a linear-quadratic problem. There, the stage cost is defined as a quadratic function and the dynamics are linear. Thus, the linear-quadratic optimal control problem is given by

\begin{equation} \begin{split} \min_{u(\cdot)\in\mathbb{U}}J^N(x_0,u(\cdot)) &= \sum_{k=0}^{N-1}\ell(x(k,x_0),u(k))+F(x(N,x_0))\\ &= \sum_{k=0}^{N-1}x(k,x_0)^T Q x(k,x_0) +u(k)^T R u(k)+ 2x(k,x_0)^TN u(k) \\ & \qquad \quad +x(N,x_0)^T P x(N,x_0)\\ \text{s.t.}\quad x(k+1,x_0)&=Ax(k,x_0)+Bu(k),\quad k = 0,\dots, N-1\\ x(0)&= x_0\\ x\in\mathbb{X},\quad & x(N,x_0)\in\mathbb{X}_0 \end{split} \end{equation}

where \(Q, R, N, P\) are weightening matrices and \(A, B\) the system matrices, each respectively of suitable dimension. Further, the constraints have to be also linear and of the form

\[Ex+ Fu \geq h.\]

Note

nMPyC supports a time dependent formulation of optimal control problem. Hence, all functions, as \(f, \ell, F\), can depend on the time instance \(j\).

Note

nMPyC supports also discounted optimal control problems. In the discrete case the objective is defined as

\[J^N(x_0,u(\cdot)):=\sum_{k=0}^{N-1}\beta^k\ell(x(k,x_0),u(k))\]

with \(\beta\in(0,1)\) the discount factor.

API Reference

system

A class used to define the system dynamics of an optimal control problem.

objective

A class used to define the objective of the optimal control problem.

constraints

Class used to define the constraints of the optimnal control problem.

model

Class that contains all the components of the optimal control problem.

result

Class used to store the simulation results of the MPC simulation.

nmpyc_array

Module for array definition and computation.

Examples

In addition to the information from the API References, the following examples are intended to provide guidance for implementing your own problems.

To show the different possibilities of the nMPyC package, we illustrate them with different examples.

Therefore, the chemical reactor is a nonlinear autonomous problem, the inverted pendulum is a linear quadratic problem, the heat pump is a nonlinear time-varying problem and the 2d investment problem is a discounted problem.

Chemical Reactor

We consider a single first-order, irreversible chemical reaction in an isothermal CSTR

\[A \to B.\]

The material balances and the system data are provided in [DAR11] and is given by the discrete time nonlinear model

\begin{equation*} \begin{split} c_{A}(k+1)&=c_{A}(k)+h\left(\frac{Q(k)}{V}(c_f^{A}-c_{A}(k))-k_r{c_{A}(k)}\right)\\ c_{B}(k+1)&=c_{B}(k)+h\left(\frac{Q(k)}{V}(c_f^{B}-c_{B}(k))+k_r{c_{B}(k)}\right), \end{split} \end{equation*}

in which \(c_A\geq 0\) and \(c_B\geq 0\) are the molar concentrations of \(A\) and \(B\) respectively, \(Q\leq 20\) (L/min) is the flow through the reactor and \(h\) is the sampling rate in minutes. The constants and their meanings are given in table below.

Reactor constants
		value	unit
feed concentration of \(A\)	\(c_f^{A}\)	1	mol/L
feed concentration of \(B\)	\(c_f^{B}\)	0	mol/L
volume of the reactor	\(V_R\)	10	L
rate constant	\(k_r\)	1.2	L/(mol min)
sampling rate	\(h\)	0.5	min
initial value	\((c_0^{A},c_0^B)\)	\((0.4, 0.2)\)

From this set of parameters we can compute the equilibrium \((c_e^{A},c_e^B,Q_e) = (\frac 1 2, \frac 1 2, 12)\) of the system.

To initialize the system dynamics a function that implements \(f(x,u)\), where \(x = (c_{A},c_{B})^T\) and \(u=Q\) has to be defined.

V = 10.
cf_A = 1.
cf_B = 0.
k_r = 1.2
h = 0.5

def f(x,u):
    y = nmpyc.array(2)
    y[0] = x[0] + h*((u[0]/V) *(cf_A - x[0]) - k_r*x[0])
    y[1] = x[1] + h*((u[0]/V) *(cf_B - x[1]) + k_r*x[1])
    return y

After that, the nMPyC system object can be set by calling

system = nmpyc.system(f, 2, 1, system_type='discrete')

In the next step, the objective is defined by using the stage cost given by

\begin{align*} \ell (c_{A}(k),c_{B}(k),Q(k))&=\frac 1 2\vert c_A(k)-\frac 1 2\vert^2+\frac 1 2 \vert c_B(k)-\frac 1 2\vert^2+\frac 1 2 \vert Q(k) -12 \vert^2\\ \end{align*}

Since we do not need terminal cost, we can initialize the objective directly using the following implementation.

def l(x,u):
    return 0.5 * (x[0]-0.5)**2 + 0.5 * (x[1]-0.5)**2 + 0.5 * (u[0]-12)**2

objective = nmpyc.objective(l)

In terms of the constraints we assume that

\[\begin{split}0 &\leq x_1(k) & < \infty & \quad & \text{for } k=0,\ldots,N \\ 0 &\leq x_2(k) & < \infty & \quad & \text{for } k=0,\ldots,N \\ 0 &\leq u(k) & \leq 20 & \quad & \text{for } k=0,\ldots,N-1.\end{split}\]

This can be realized in the code as follows:

constraints = nmpyc.constraints()
lbx = nmpyc.zeros(2)
ubu = nmpyc.ones(1)*20
lbu = nmpyc.zeros(1)
constraints.add_bound('lower','state', lbx)
constraints.add_bound('lower','control', lbu)
constraints.add_bound('upper','control', ubu)

Moreover, we consider the equilibrium \((c_e^{A},c_e^B,Q_e)\) as th terminal condition for our optimal control problem, which is implemented as

xeq = nmpyc.array([0.5,0.5])
def he(x):
    return x - xeq
constraints.add_constr('terminal_eq', he)

After all components of the optimal control problem have been implemented, we can now combine them into a model and start the MPC loop. For this Purpose, we define

\[x(0) = (0.4,0.2)^T\]

and set \(N=15\), \(K=100\).

model = nmpyc.model(objective,system,constraints)
x0 = nmpyc.array([0.4,0.2])
res = model.mpc(x0,15,100)

Following the simulation we can visualize the results by calling

res.plot()

which generates the plot bellow.

Inverted Pendulum

We consider the mechanical model of an inverted rigid pendulum mounted on a carriage, see [Grune21], [GruneP17].

By means of physical laws an “exact” differential equation model can be derived. However, since in our case we like to obtain a linear quadratic problem, we linearize the differential equation at the origin. Thus, we obtain the system dynamics defined by

\[\begin{split}\dot{x}(t) = \left(\begin{array}{cccc} 0 & 1 & 0 & 0 \\ g & -k & 0 & 0 \\ 0 & 0 & 0 & 1 \\ 0 & 0 & 0 & 0 \end{array}\right) x(t) + \left(\begin{array}{c} 0 \\ 1 \\ 0 \\ 1 \end{array}\right) u(t).\end{split}\]

Here, the state vector \(x \in \mathbb{R}^4\) consists of 4 components. \(x_1\) corresponds to the angle \(\psi\) of the pendulum, which increases counterclockwise, where \(x_1 = 0\) corresponds to the upright pendulum. \(x_2\) is the angular velocity, \(x_3\) the position of the carriage and \(x_4\) its velocity. The control \(u\) is the acceleration of the carriage. The constant \(k=0.1\) describes the friction of the pendulum and the constant \(g \approx 9.81 m/s^2\) is the acceleration due to gravity.

Since the system dynamics are linear, we can initialize them using the LQP method.

g = 9.81
k = 0.1
A = nmpyc.array([[0, 1, 0, 0],
                 [g, -k, 0, 0],
                 [0, 0, 0, 1],
                 [0, 0, 0, 0]])
B = nmpyc.array([0, 1, 0, 1])
system = nmpyc.system.LQP(A, B, 4, 1, 'continuous',
                          sampling_rate=0.1, method='rk4')

Note that we have to use one of the fixed step methods as euler, heun or rk4 as integration method if we like to exploit the linear quadratic structure of the problem in the optimization.

In the next step, we have to define the objective of the optimal control problem. In doing so, we assume the stage cost

\[\ell(x,u) = 2x^Tx + 4u^Tu.\]

Since we assume no terminal cost, we can implement the objective as shown in the following code snippet.

Q = 2*nmpyc.eye(4)
R = 4*nmpyc.eye(1)
objective = nmpyc.objective.LQP(Q, R)

Again, we use the LQP method to exploit the linear quadratic structure of the problem later.

In terms of the constraints we assume the state constraints

\[-9 \leq x_i(t) \leq 5\]

for \(i=1,\ldots,4\) and the control constraint

\[-20 \leq u(t) \leq 6 \quad\]

This can be realized in the code as

constraints = nmpyc.constraints()
lbx = nmpyc.zeros(4)*(-9)
ubx = nmpyc.ones(4)*5
constraints.add_bound('lower','state', lbx)
constraints.add_bound('upper','state', ubx)
constraints.add_bound('lower', 'control', nmpyc.array([-20]))
constraints.add_bound('upper', 'control', nmpyc.array([6]))

After all components of the optimal control problem have been implemented, we can now combine them into a model and start the MPC loop. For this Purpose, we define the inital value

\[x(0) = (1,1,1,1)^T\]

and set \(N=20\), \(K=100\).

model = nmpyc.model(objective,system,constraints)
x0 = nmpyc.array([1, 1, 1, 1])
res = model.mpc(x0,20,100)

Since the problem is linear-quadratic, the program automatically takes advantage of this fact and uses the appropriate solver osqp. To change this and use for example the SciPy solver SLSQP, we can use the set_options method before calling model.mpc().

model.opti.set_options(dict(solver='SLSQP'))

Note that changing the optimizer usually does not have any advantage and is therefore not necessarily recommended. At this point we only like to demomnstrate the use of this function.

Following the simulation we can visualize the open and closed loop results by calling

res.plot() # plot closed loop results
res.plot('state', show_ol=True) # plot open loop states

which generates the plots bellow.

Heat Pump

This example describes a home heating system that involves the optimal control of a small heat pump coupled to a floor heating system. The corresponding dynamic model is introduced in [LHDI10] and is given by

\begin{align} \dot{x_1} &= \dfrac{-k_{WR}}{\rho_W c_W V_H}x_1 + \dfrac{k_{WR}}{\rho_W c_W V_H}x_2 + \dfrac{1}{\rho_W c_W V_H}u\\ \dot{x_2}&= \dfrac{k_{WR}}{k_G \tau_G}x_1 -\dfrac{k_{WR}+k_G}{k_G\tau_G}x_2 +\dfrac{1}{\tau_G} T_{\text{amb}}, \end{align}

where \(x_1\) denotes the temperature of the water returning from the heating, \(x_2\) denotes the room temperature and \(u\) is the heat supplied from the heat pump to the floor. Further, the ambient temperature

\[T_\text{amb}(t) = 2.5 + 7.5 \sin\left(\frac{2\pi t}{t_f}-\frac \pi 2\right)\]

describes a sinusoidal disturbance from the outside temperature where \(t_f = 24\). The remaining constants are summarized in the table below.

Reactor constants
		value	unit
density of the water	\(\rho_W\)	997	\(kg/m^3\)
specific heat capacity of water	\(c_W\)	4.1851	\(J/kgK\)
volume of the water	\(V_H\)	7.4	\(m^3\)
thermal conductivity between water and the room	\(k_{WR}\)	510	\(W/K\)
thermal conductivity between the room and the environment	\(k_G\)	125	\(W/K\)
thermal time constant of the room	\(\tau_G\)	260	\(s\)

First, we have to implement the outside temperature in the code to define our system dynamics.

t_f = 24
def T_amb(t):
 return 2.5 + 7.5*nmpyc.sin((2*nmpyc.pi*t)/t_f - (nmpyc.pi/2))

After that, we can define the right hand side of the system by

rho_W = 997
c_W = 4.1851
V_H = 7.4
k_WR = 510
k_G = 125
thau_G = 260
def f(t,x,u):
   y = nmpyc.array(2)
   y[0] = (-k_WR/(rho_W*c_W*V_H)*x[0]
            + k_WR/(rho_W*c_W*V_H)*x[1]
            + 1/(rho_W*c_W*V_H)*u[0])
   y[1] = (k_WR/(k_G*thau_G)*x[0]
            - (k_WR + k_G)/(k_G*thau_G)*x[1]
            + (1/thau_G)*T_amb(t))
   return y

And finally initialize the system by

system = nmpyc.system(f, 2, 1, 'continuous', sampling_rate=0.5, method='euler')

In the heating system the conflict between energy and thermal comfort arises. Thus, the stage cost reads

\begin{align*} \ell(x,u)&=\frac u P_{\max} + (x_2-T_\text{ref})^2, \end{align*}

where \(P_{\max} = 15000 (W)\) is the maximal power of the heating pump and \(T_\text{ref} = 22^{\circ} C\) is the desired temperature of the room. The reference temperature \(T_\text{ref}\) can be selected differently – depending on the individual thermal comfort.

According to this we can initialize our objective by

P_max = 15000
T_ref = 22
def l(x,u):
   return (u[0]/P_max) + (x[1]-T_ref)**2

and implement the control constraint

\[0 \leq u(t) \leq P_{max}\]

as

constraints = nmpyc.constraints()
constraints.add_bound('lower', 'control', nmpyc.array([0]))
constraints.add_bound('upper', 'control', nmpyc.array([P_max]))

After all components of the optimal control problem have been implemented, we can now combine them into a model and start the MPC loop. For this purpose, we define

\[x(0) = (22, 19.5)^T\]

and set \(N=30\) and \(K=500\).

model = nmpyc.model(objective,system,constraints)
x0 = nmpyc.array([22., 19.5])
res = model.nmpyc(x0,N,K)

Following the simulation we can visualize the results by calling

res.plot()

which generates the plot bellow.

2d Investment Problem

To examplify a discounted problem we consider a 2d variant of an investment problem, originally introduced in [HKHF03] and furhter explored in [GruneSS15].

The system dynamics for this problem are given by

\[\begin{split}\dot{x}_1(t) &= x_2(t) - \sigma x_1(t) \\ \dot{x}_2(t) &= u(t)\end{split}\]

where we set \(\sigma = 0.25\) for our calculations.

To implement the dynamics we have to initialize a function that implements the right hand side of the dynamics.

sigma = 0.25
def f(x,u):
    y = nmpyc.array(2)
    y[0] = x[1]-sigma *x[0]
    y[1] = u
    return y

After that, the nMPyC system object can be set by calling

system = nmpyc.system(f, 2, 1, 'continuous', sampling_rate=0.2, method='heun')

To model the payoff of the investment problem we assume th stage cost

\[\ell(x,u) = -(R(x_1) - c(x_2) - v(u))\]

where \(R(x_1) = k_1 \sqrt{x_1} - x_1/(1+k_2 x_1^4)\) is a revenue function of the firm with a convex segment due to increasing returns. \(c(x_2) = c_1 x_2 + c_2 x_2^2/2\) denotes adjustment costs of investment and \(v(u) = \alpha u^2/2\) represents adjustment costs of the change of investment. The convex segment in the payoff function just mentioned is likely to generate two domains of attraction. Additionally we choose \(k_1=2\), \(k_2=0.0117\), \(c_1=0.75\), \(c_2=2.5\) and \(\alpha=12\) for our computations.

With the nMPyC package the implemnetiation of the objective corresponding to this costs can be done as follws.

def l(x,u):
    R = k1*x[0]**(1/2)-x[0]/(1+k2*x[0]**4)
    c = c1*x[1]+(c2*x[1]**2)/2
    v = (alpha*u[0]**2)/2
    return -(R - c - v)

objective = nmpyc.objective(l)

Since this problem is unconstrained we can now initialize our model by

model = nmpyc.model(objective,system)

For our simulation we assume set the discount factor to

\[\beta = e^{-\delta h}\]

where \(h=0.2\) is our samplimng rate and \(\delta=0.04\) is the continuous discount rate.

It can now be shown that this problem has two domains of attraction, one at roughly \(x^* = (0.5, 0.2)\) and the other roughly at \(x^* = (4.2, 1.1)\). Now we choose ifferent initial values from both domains of attraction to test, if we can replicate the two domains of attraction for a finite decision horizon by using nonlinear model predictive control. For this purpose we set the MPC horizon \(N=50\) and the number of MPC iterations to \(K=500\).

This leads to the following code for running the closed loop simulation for the discounted problem.

discount = nmpyc.exp(-0.04*0.2)
N = 50
K = 500

x0 = nmpyc.array([3.0,0.75])
res1 = nmpyc.mpc(x0,N,K,discount)

x0 = nmpyc.array([5.0,1.75])
res2 = nmpyc.mpc(x0,N,K,discount)

Looking at the phase portraits of the two simulations, we can confirm that we really converge against the two different equilibria with the closed loop trajectory. The phase portraits of our simulations can be plotted with the nMPyC package by calling

res1.plot('phase', phase1='x_1', phase2='x_2', show_ol=True)
res2.plot('phase', phase1='x_1', phase2='x_2', show_ol=True)

The option show_ol=True will also plot the pahase portraits of the open loop simulations of each iteration, which leads the output below.

Templates

In addition to the examples, we also provide templates to facilitate the implementation.

To take advantage of the different structures of the problems, we have implemented templates for the following problem types.

Time-variant Problem

# Import nMPyc package
import nmpyc

# Define system parameters
nx = .. # dimension of state
nu = .. # dimension of control
system_type = .. # system type: continuous or discrete
sampling_rate = 1. # sampling rate h (optional)
t0 = 0. # initial time (optional)
method = 'cvodes' # integrator (optinal)

# Define MPC parameters
N = .. # MPC horizon
K = .. # MPC itertaions
x0 = .. # initial value
discount = 1. # dicount factor (optional)

# Define right hand side of the system dynamics
def f(t, x, u):
   y = nmpyc.array(nx)
   ..
   return y

# Initialize system dynamics
system = nmpyc.system(f, nx, nu, system_type, sampling_rate, t0, method)

# Define stage cost
def l(t, x, u):
   return ..

# Define terminal cost (optional)
def F(t, x):
   return ..

# Initialize objective
objective = nmpyc.objective(l, F)

# Define constraints
constraints = nmpyc.constraints()

# Add bounds (optional)
lbx = .. # lower bound for states
constraints.add_bound('lower', 'state', lbx)
ubx = .. # upper bound for states
constraints.add_bound('upper', 'state', ubx)
lbu = .. # lower bound for control
constraints.add_bound('lower', 'control', lbu)
ubu = .. # upper bound for control
constraints.add_bound('upper', 'control', ubu)
lbend = .. # lower bound for terminal state
constraints.add_bound('lower', 'terminal', lbend)
ubend = .. # upper bound for terminal state
constraints.add_bound('upper', 'terminal', ubend)

# Add equality constraints (h(t,x,u)=0, optional)
len_eqconstr = .. # number of equality constraints
def h(t, x, u):
   c_eq = nmpyc.array(len_eqconstr)
   ..
   return c_eq

constraints.add_constr('eq', h)

# Add inequality constraints (g(t,x,u)>=0, optional)
len_ineqconstr = .. # number of inequality constraints
def g(t, x, u):
   c_ineq = nmpyc.array(len_ineqconstr)
   ..
   return c_ineq

constraints.add_constr('ineq', g)

# Add terminal equality constraints (H(t,x)=0, optional)
len_terminaleq = .. # number of terminal equality constraints
def H(t, x):
   cend_eq = nmpyc.array(len_terminaleq)
   ..
   return cend_eq

constraints.add_constr('terminal_eq', H)

# Add terminal equality constraints (G(t,x)>=0, optional)
len_terminalineq = .. # number of terminal equality constraints
def G(t, x):
   cend_ineq = nmpyc.array(len_terminalineq)
   ..
   return cend_ineq

constraints.add_constr('terminal_ineq', G)

# Initialize model
model = nmpyc.model(objective, system, constraints)

# Start MPC loop
res = model.mpc(x0, N, K, discount)

# Plot results
res.plot()

Autonomous Problem

# Import nMPyc package
import nmpyc

# Define system parameters
nx = .. # dimension of state
nu = .. # dimension of control
system_type = .. # system type: contiunous or discrete
sampling_rate = 1. # sampling rate h (optional)
method = 'cvodes' # integrator (optinal)

# Define MPC parameters
N = .. # MPC horizon
K = .. # MPC itertaions
x0 = .. # initial value
discount = 1. # dicount factor (optional)

# Define right hand side of the system dynamics
def f(x, u):
   y = nmpyc.array(nx)
   ..
   return y

# Initialize system dynamics
system = nmpyc.system(f, nx, nu, system_type, sampling_rate, method=method)

# Define stage cost
def l(x, u):
   return ..

# Define terminal cost (optional)
def F(x):
   return ..

# Initialize objective
objective = nmpyc.objective(l, F)

# Define constraints
constraints = nmpyc.constraints()

# Add bounds (optional)
lbx = .. # lower bound for states
constraints.add_bound('lower', 'state', lbx)
ubx = .. # upper bound for states
constraints.add_bound('upper', 'state', ubx)
lbu = .. # lower bound for control
constraints.add_bound('lower', 'control', lbu)
ubu = .. # upper bound for control
constraints.add_bound('upper', 'control', ubu)
lbend = .. # lower bound for terminal state
constraints.add_bound('lower', 'terminal', lbend)
ubend = .. # upper bound for terminal state
constraints.add_bound('upper', 'terminal', ubend)

# Add equality constraints (h(x,u)=0, optional)
len_eqconstr = .. # number of equality constraints
def h(x, u):
   c_eq = nmpyc.array(len_eqconstr)
   ..
   return c_eq

constraints.add_constr('eq', h)

# Add inequality constraints (g(x,u)>=0, optional)
len_ineqconstr = .. # number of inequality constraints
def g(x, u):
   c_ineq = nmpyc.array(len_ineqconstr)
   ..
   return c_ineq

constraints.add_constr('ineq', g)

# Add terminal equality constraints (H(x)=0, optional)
len_terminaleq = .. # number of terminal equality constraints
def H(x):
   cend_eq = nmpyc.array(len_terminaleq)
   ..
   return cend_eq

constraints.add_constr('terminal_eq', H)

# Add terminal equality constraints (G(x)>=0, optional)
len_terminalineq = .. # number of terminal equality constraints
def G(x):
   cend_ineq = nmpyc.array(len_terminalineq)
   ..
   return cend_ineq

constraints.add_constr('terminal_ineq', G)

# Initialize model
model = nmpyc.model(objective, system, constraints)

# Start MPC loop
res = model.mpc(x0, N, K, discount)

# Plot results
res.plot()

Linear Quadratic Problem

# Import nMPyc package
import nmpyc

# Define system parameters
nx = .. # dimension of state
nu = .. # dimension of control
system_type = .. # system type: continuous or discrete
sampling_rate = 1. # sampling rate h (optional)
method = 'cvodes' # integrator (optinal)

# Define MPC parameters
N = .. # MPC horizon
K = .. # MPC itertaions
x0 = .. # initial value
discount = 1. # dicount factor (optional)

# Define linear right hand side of the system dynamics f(x,u) = Ax + Bu
A = ..
B = ..

# Initialize system dynamics
system = nmpyc.system.LQP(A, B, nx, nu, system_type, sampling_rate, method=method)

# Define quadratic stage cost l(x,u) = x^TQx + u^TRu + 2*x^TNx
Q = ..
R = ..
N = nmpyc.zeros((nx,nu)) # optional

# Define terminal cost x^TPx
P = nmpyc.zeros((nx,nx)) # optional

# Initialize objective
objective = nmpyc.objective.LQP(Q, R, N, P)

# Define constraints
constraints = nmpyc.constraints()

# Add bounds (optional)
lbx = .. # lower bound for states
constraints.add_bound('lower', 'state', lbx)
ubx = .. # upper bound for states
constraints.add_bound('upper', 'state', ubx)
lbu = .. # lower bound for control
constraints.add_bound('lower', 'control', lbu)
ubu = .. # upper bound for control
constraints.add_bound('upper', 'control', ubu)
lbend = .. # lower bound for terminal state
constraints.add_bound('lower', 'terminal', lbend)
ubend = .. # upper bound for terminal state
constraints.add_bound('upper', 'terminal', ubend)

# Add equality constraints (Ex + Fu = b, optional)
E_eq = ..
F_eq = ..
b_eq = ..
constraints.add_constr('eq', E_eq, F_eq, b_eq)

# Add equality constraints (Ex + Fu >= b, optional)
E_ineq = ..
F_ineq = ..
b_ineq = ..
constraints.add_constr('ineq', E_ineq, F_ineq, b_ineq)

# Add terminal equality constraints (Hx = 0, optional)
H_eq = ..
constraints.add_constr('terminal_eq', H_eq)

# Add terminal equality constraints (Hx >= 0, optional)
H_ineq = ..
constraints.add_constr('terminal_ineq', H_ineq)

# Initialize model
model = nmpyc.model(objective, system, constraints)

# Start MPC loop
res = model.mpc(x0, N, K, discount)

# Plot results
res.plot()

Note

Any problem, whether nonlinear, linear, autonomous, or time-varying, can be initialized as a nonlinear time-varying optimal control problem. Therefore, you can always fall back on such an implementation. However, if you know the structure of your problem and this is to be exploited by the program in order to possibly speed up the simulation, it is necessary to initialize the problem as such.

FAQ

Why use the nmpyc_array module?

The idea of the nmpyc.nmpyc_array module is to provide a simple syntax for the input, which is similar to the one of NumPy. At the same time we ensure a switching between symbolic calculation with CasADi and completely numeric calculations.

A completely numerical calculation is advantageous, for example, if non-differentiable functions have to be evaluated at critical points, e.g. the norm at the origin. Here the algorithmic differentiation of CasADi can lead to problems.

Therefore this module is built in a way that the array class can automatically switch between CasADi and NumPy objects. In addition, the individual functions are built in such a way that they recognize the type of the input and call the appropriate function from NumPy or CasADi accordingly.

What to do if a function is not defined in the nmpyc_array module?

We have tried to implement the most common functions. Nevertheless, it can happen that a certain function that you need is missing.

If this is the case, there is the possibility to implement an own overload of this function. A good orientation for this is the already programmed functions in the nmpyc.nmpyc_array module.

Another possibility in such cases is to use the call x.A to access the CasADi or NumPy array in which the entries of x are stored. Afterwards the appropriate necessary computations can be accomplished with the help of NumPy or CasADi functions. Note, however, that in this way if applicable no smooth change between numeric and symbolic calculation is possible.

Which solver should be used?

In the most cases, the automatic selection of the solver by the program is recommended. In this way, if possible, the linear quadratic structure of a problem is exploited or at least algorithmic differentiation is still exploited to perform an advantageous optimization.

However, as already mentioned, this algorithmic differentiation can also lead to problems in some cases. For example, if a non-differentiable function must be evaluated at critical points, e.g. the norm at the origin. In such cases, a numerical calculation should be used for the optimization and a SciPy solver, such as SLSQP, should be selected.

Which discretization method should be used?

In our numerical simulations we have experienced that mostly a fixed step integration method like euler is sufficient to guarantee the necessary accuracy during the simulation. The advantage of these methods is that with them the largest speed up among the available integrators can be achieved.

However, if it is necessary to achieve higher integration accuracy by an adaptive integration method, one of the CasADi integrators, e.g. cvodes, should always be chosen if possible.

The SciPy integrators should only be considered as a kind of backup in case the other methods fail, since they lead to an above-average lag of time during the simulation in our implementation.

What to do if I a LaTeX Error occurs while plotting?

In our experience such errors occur mainly on MacOS if Spyder is used for programming, which in turn is opened via the Anaconda Navigator. In this case it is sufficient to open spyder directly and not to take the detour via the Anaconda Navigator to solve the problem.

However, if this procedure does not solve the problem or the problem has another cause, it is also possible to disable the LaTeX labeling of the plots by setting the option usetex=False. For more details see nmpyc.result.result.plot().

How to Cite

If you use nMPyC for published work please cite it as

@misc{nmpyc,
       author = {Jonas Schie{\ss}l and Lisa Kr{\"u}gel},
       title ={{nMPyC} - A Python library for solving optimal control problems via MPC},
       howpublished = {\url{http://nmpyc.readthedocs.io/}},
       year = {2022}
}

Please remember to properly cite other software that you might be using too if you use (e.g. CasADi, IPOPT, …).

For any specific algorithm, also consider citing the original author’s paper.

References

[DAR11]

Moritz Diehl, Rishi Amrit, and James B. Rawlings. A lyapunov function for economic optimizing model predictive control. IEEE Transactions on Automatic Control, 56(3):703–707, mar 2011. doi:10.1109/TAC.2010.2101291.

[Grune21]

Lars Grüne. Mathematical control theory. 2021. Lecture Notes. URL: https://num.math.uni-bayreuth.de/de/team/lars-gruene/skripten/kontrolltheorie/kt_2021_en.pdf.

[GruneP17]

Lars Grüne and Jürgen Pannek. Nonlinear Model Predictive Control : Theory and Algorithms. 2nd Edition. Communications and Control Engineering. Springer, Cham, Switzerland, 2017. URL: https://eref.uni-bayreuth.de/35127/.

[GruneSS15]

Lars Grüne, Willi Semmler, and Marleen Stieler. Using nonlinear model predictive control for dynamic decision problems in economics. Journal of Economic Dynamics and Control, 60:112–133, 2015. URL: https://eref.uni-bayreuth.de/20841/.

[HKHF03]

Josef L. Haunschmied, Peter M. Kort, Richard F. Hartl, and Gustav Feichtinger. A DNS-curve in a two-state capital accumulation model: a numerical analysis. Journal of Economic Dynamics and Control, 27(4):701–716, feb 2003. doi:10.1016/S0165-1889(01)00070-7.

[LHDI10]

Filip Logist, Boris Houska, Moritz Diehl, and Jan Van Impe. Fast pareto set generation for nonlinear optimal control problems with multiple objectives. Structural and Multidisciplinary Optimization, 42(4):591–603, may 2010. doi:10.1007/s00158-010-0506-x.

Welcome to nMPyC’s documentation!

Contents

Installation

Requirements

Installation using PIP

Installation by Source

Getting Started

Import nMPyC

Creating the System Dynamics

Creating the Objective

Creating the Constraints

Running Simulations

Advanced topics

Basics of model predictive control

Optimal control problems

The basic MPC algorithm

Notes and extensions

Further reading

API Reference

Examples

Chemical Reactor

Inverted Pendulum

Heat Pump

2d Investment Problem

Templates

Time-variant Problem

Autonomous Problem

Linear Quadratic Problem

FAQ

Why use the nmpyc_array module?

What to do if a function is not defined in the nmpyc_array module?

Which solver should be used?

Which discretization method should be used?

What to do if I a LaTeX Error occurs while plotting?

How to Cite

References