hilo-mpc / hilo-mpc

HILO-MPC is a Python toolbox for easy, flexible and fast development of machine-learning-supported optimal control and estimation problems

HILO-MPC

HILO-MPC is a Python toolbox for easy, flexible and fast realization of machine-learning-supported optimal control, and estimation problems developed mainly at the Control and Cyber-Physical Systems Laboratory, TU Darmstadt, and the Laboratory for Systems Theory and Control, Otto von Guericke University. It can be used for model predictive control, moving horizon estimation, Kalman filters, solving optimal control problems and has interfaces to embedded model predictive control tools.

HILO-MPC can interface directly to TensorFlow and PyTorch to create machine learning models and the CasADi framework to efficiently build control and estimation problems. The machine learning models can be used (almost) everywhere in the setup of these problems.

Currently the following machine learning models are supported:

• Feedforward neural networks
• Gaussian processes

The following machine learning models are currently under development:

• Bayesian neural network
• Recurrent neural network

At the moment the following MPC and optimal control problems can be solved:

• Reference tracking nonlinear MPC
• Trajectory tracking nonlinear MPC
• Path following nonlinear MPC
• Economic nonlinear MPC
• Linear MPC
• Traditional optimal control problems

All the nonlinear MPCs support soft constraints, time-variant systems, time-varying parameters and can be used to solve minimum-time problems. They work for continuous-time and discrete-time models, in DAE or ODE form. Linear MPC is currently limited towards discrete-time models.

A rich set of examples is available, spanning:

• NMPC for bioreactors using hybrid first principle and learned models
• Trajectory tracking and path following model predictive control with learning and obstacle avoidance
• Output feedback MPC of a continuous stirred tank reactor with a Gaussian process prediction model
• Learning NMPC control using a neural network
• Simple LQR, PID
• Moving horizon estimation, extended Kalman filter, unscented Kalman filter, and particle filter for a continuous stirred tank reactor

Installation

You can use pip to install HILO-MPC as follows

pip install hilo-mpc

Additional Packages

If you want to make use of the complete functionality of the toolbox, you may want to install one of the following packages

Package	Version	Usage
TensorFlow	>=2.3.0, <2.8.0	Training of neural networks
PyTorch	>=1.2.0	Training of neural networks
scikit-learn	>=0.19.2	Normalizing of training data
Bokeh	>=2.3.0, <3.0.0	Plotting
Matplotlib	>=3.0.0	Plotting
pandas	>=1.0.0 <=1.5.1	Providing data for training of neural networks

Documentation

A preliminary documentation can be found here. Note that this documentation is not complete and will be updated over time.

Citing HILO-MPC

If you use HILO-MPC for your research, please cite the following publication:

• J. Pohlodek, B. Morabito, C. Schlauch, P. Zometa, R. Findeisen. Flexible development and evaluation of machine-learning-supported optimal control and estimation methods via HILO-MPC. arXiv. 2022.

@misc{pohlodek2022hilompc,
    title = {Flexible development and evaluation of machine-learning-supported optimal control and estimation methods via {HILO-MPC}},
    author = {Pohlodek, Johannes and Morabito, Bruno and Schlauch, Christian and Zometa, Pablo and Findeisen, Rolf},
    publisher = {arXiv},
    year = {2022},
    doi = {10.48550/ARXIV.2203.13671}
}

hilo-mpc's Issues

Allow the user to easily obtain the values of stage and terminal cost at every iteration

          Hi @aabsz,

unfortunately there's no direct way how to get that value, but what you can do is, since you know how these terms looks like, you can create a function that reproduce them and evaluate these functions at every time step

Originally posted by @brunomorampc in #41 (comment)

Improve plot_iterations

Bugs:

1. not working if the number of iterations is large

If the number of iterations it largen than 256, the palette inferno will return an error.

Solution

Use same color OR use a larger palette OR loop over the palette again.

Improvements

1. Add legend to section shouwing the different kind of constraints

Ability to add both inputs and input_change cost to quad cost

I've noticed that when you use quad_stage_cost it is not possible to use both add_inputs and add_inputs_change for the same set of inputs which is commonly required in many MPC problems. It results in the error:

TypeError("Two different varying trajectory for the same states are not allowed.")

I got around that by modifying _check_options in modeling.py like so:

if check_if_has_duplicates(self.name_varying_trajectories) and references is not None:
            raise TypeError("Two different varying trajectory for the same states are not allowed.")

I enabled IPOPT output and tested on the problem I am trying to solve, the objective values do change when I use add_inputs and add or remove add_inputs_change, and when I run the controller for a long trajectory I can see a difference in the solutions, the addition of add_inputs_change does smooth the control actions.

I don't know if this change breaks anything else in the library or if it is the right solution, so I thought I'd bring it up here as it would be nice to enable both add_inputs and add_inputs_change. What are your thoughts on that? Am I breaking something else by modifying that line?

On a separate note, I have noticed a huge speed up when I pass "expand": True to _nlp_opts, it would be useful to have that as an optional argument on NMPC.setup().

p.s. This library is awesome!

Problem with defining variable parameters for terminal cost function

Hi,

I am trying to define a terminal cost function with variable parameters. I know how to define variable parameters and I have already done it when defining a stage cost function with variable parameters and it is working well (using stage_cost.cost).

However, when defining the terminal cost using the following command, I get the error that the parameters have not been declared:

self.nmpc.quad_terminal_cost.add_states(names=['x', 'y', 'z'], weights=[self.model.p[1] , self.model.p[2], self.model.p[3]],ref=[self.x_dess,self.y_dess, self.z_dess])

This is the error that I get:

Initialization failed since variables [stage_cost_weight_xstage_cost_weight_ystage_cost_weight_z] are free. These symbols occur in the output expressions but you forgot to declare these as inputs. Set option 'allow_free' to allow free variables.

Is there any option that I can define the parameter varying terminal cost function?

Can hilo-mpc train the ODE model parameters according to the actual value of the ODE model parameters without knowing

I studied example: nmpc_hybrid_bio.ipynb. Then I found that the ANN model fits 'mu', 'Rs', 'Rfp' according to 'S', 'I' of the original data. But in general, the specific reaction rate of the fermentation process is unknown. How should I use hilo-mpc to train the kinetic parameters based on the existing data (including t, X, S, P, I, V, Feed, etc.) model, and then combine the two to build a Hybrid model.

Optimization fails because x0 is None.

I'm trying to implement a NMPC for trajectory tracking with actual reference values(coordinates) instead of a function.

set_plot_backend('bokeh')

# Define the model
model = Model(name='MPC trajectory', plot_backend='bokeh', discrete=True)

#States
X = model.set_dynamical_states(['x', 'y', 'theta', 'v'])
model.set_measurements(['y_x', 'y_y', 'y_theta', 'y_v'])
model.set_measurement_equations([X[0], X[1], X[2], X[3]])

#Unwrap the states
x = X[0]
y = X[1]
theta = X[2]
v = X[3]
#Control Inputs
U = model.set_inputs(['a', 'delta'])

#Unwrap the inputs
a = U[0]
delta = U[1]

#Kinematic bicycle model
lr = 1.5
lf = 1.5
beta = ca.arctan(lr / (lr + lf) * ca.tan(delta))
w = v * ca.tan(delta) * ca.cos(beta) / (lr + lf)
dx = v * ca.cos(theta + beta)
dy = v * ca.sin(theta + beta)
dtheta = w
dv = a
model.set_dynamical_equations([dx, dy, dtheta, dv])

#Initial conditions
x0 = [-6.446169853210449, -79.05502319335938, 0, 0]
u0 = [0,0]

model.setup(options={'objective_function': 'discrete'})
nmpc = NMPC(model)

#List of some reference values
x_ref = [.....]
y_ref = [.....]

nmpc.quad_stage_cost.add_states(names=['x', 'y'], weights=[10., 10.], trajectory_tracking=True)
nmpc.quad_terminal_cost.add_states(names=['x', 'y'], weights=[10.,10.], trajectory_tracking=True)
nmpc.quad_stage_cost.add_states(names=['v'], ref=2, weights=10)
nmpc.quad_terminal_cost.add_states(names=['v'], ref=2, weights=10)

nmpc.horizon = 10
nmpc.set_initial_guess(x_guess=x0, u_guess=u0)
nmpc.setup(options={'objective_function': 'discrete', 'print_level': 0})

ref_sc_dict = {'x': x_ref, 'y': y_ref}
ref_tc_dict = {'x': x_ref, 'y': y_ref}

model.set_initial_conditions(x0=x0)
for step in range(100):                                                     
     u = nmpc.optimize(xt, ref_sc=ref_sc_dict, ref_tc=ref_tc_dict)
     model.simulate(u=u)
     x0 = model.solution['xf']

I keep running into this error that indicates model.solution['xf'] has no values.

Traceback (most recent call last):
  File ".trial.py", line 139, in <module>
    u = nmpc.optimize(x0, ref_sc=ref_sc_dict, ref_tc=ref_tc_dict)
  File "hilo_mpc\modules\controller\mpc.py", line 777, in optimize
    x0 = check_and_wrap_to_DM(x0)
  File "hilo_mpc\util\util.py", line 236, in check_and_wrap_to_DM
    raise TypeError(f"Type {type(arg)} not supported. Must be list,float,ca.DM or numpy array ")
TypeError: Type <class 'NoneType'> not supported. Must be list,float,ca.DM or numpy array

Is this a bug or am I making any mistakes? Please help @brunomorampc and thanks for this amazing package!

How to use tht function nmpc.stage_cost

My objective function is non quadratic, so I want to use the function "nmpc.stage_cost". But I didn't find the relevant guide.
Does anyone know how to use this function?
Or does anyone know how to establish a non quadratic objective function by other methods?

Implement infinite-horizon LQR

Add infinite-horizon LQR to the toolbox

Defining a cost function with some terms not including model states ans inputs

Hi,

I'm currently working on a problem where my cost function has some terms that aren't directly tied to the model inputs and states. Unfortunately, I can't define these terms using the standard state-space equation format (x_dot = Ax + Bu) to include them in the model.

Specifically, I'm tackling the issue of collision avoidance for a robot. My cost function consists of two components. The first component is related to achieving a desired goal, which can be defined using the model states. However, the second component is all about avoiding collisions with obstacles. To accomplish this, I need to maximize the distance between the robot and the obstacle, adding a term proportional to (1/dist) to the cost function. The challenge here is that this distance isn't defined as a model variable with a dynamical equation, causing the cost function to ignore it.

I'm open to any suggestions you might have on how to effectively incorporate these non-standard terms into my cost function and make them influence the optimization process. Your insights would be greatly appreciated.

Getting the value of stage_cost and terminal_cost

Hi,

Would you please let me know how I can get the value of stage_cost and terminal_cost at each time step of the simulation?

MHE - Add posibility to estimate only specific model parameters

At the moment all the parameters of the model are estimated and the only way to fix some parameters (e.g. because they're known) is to put the lower and upper bound of the MHE optimization to the same value.

Add a method where the parameters can be set as known and hence they do not need to be estimated.

Prediction and control horizon

When I only use nmpc.horizon = 10, everything is right.
But when I want to set different control and prediction horizon, i.e.
nmpc.prediction_horizon=10 // nmpc.control_horizon=4
there exists an error:
`Traceback (most recent call last):

File ~/plexe-pyapi/examples/untitled0.py:42 in
nmpc.setup(options={'print_level':0})

File ~/anaconda3/lib/python3.8/site-packages/hilo_mpc/modules/controller/mpc.py:1615 in setup
u_ii = u[ii, 0]

IndexError: index 4 is out of bounds for axis 0 with size 4`

What's the problem, please

Add Bayesian NN

Add Bayesian NNs to the toolbox

Allow to add box constraints by name

At the moment the user must provide all upper bounds and lower bounds of all variables, only if just one variable is constrained.
Allow to specify the bounds of specific variables by name. Something like

mpc.set_box_constraints(names=['vx', 'vy'], ub = [10, 7] , lb = [-10, -7])

Side notes

This would decrease the computational burthen of approaches e.g. SMPC

Solve a deprecated warning in the integrator function for CasADi 3.6.3

The following deprecation warning occour with the version 3.6.3 of CasADi.

The same functionality is provided by providing additional input arguments to the 'integrator' function, in particular:
 * Call integrator(..., t0, tf, options) for a single output time, or
 * Call integrator(..., t0, grid, options) for multiple grid points.
The legacy 'output_t0' option can be emulated by including or excluding 't0' in 'grid'.
Backwards compatibility is provided in this release only.") [.../casadi/core/integrator.cpp:515]```

In [this line](https://github.com/hilo-mpc/hilo-mpc/blob/master/hilo_mpc/modules/base.py#L1713). 
It seems that the t0 and tf and grid needs to be passed as argument and not anymore in the options. 

We need to decide if we want to abandon the other version of casadi and update to the new one or allow compatibility with older versions

Update documentation

See issue #20

Particle filter initialization error

I'm trying to initialize the particle filter in the following way:

x0 = [0, 0, 0]
dt = 0.001
model=initialize(x0,dt)
observer= PF(model, roughening=True, plot_backend='bokeh')
observer.setup()

my model is Lorenz attractor initialized in the following way:

def initialize(x0,dt):

    model = Model(plot_backend='bokeh')

    states = model.set_dynamical_states(['x1', 'x2', 'x3'])
    inputs = model.set_inputs(['u1','u2','u3'])
    model.set_measurements(['o1', 'o2', 'o3'])
    model.set_measurement_equations(states)

    # Unwrap states
    x1 = states[0]
    x2 = states[1]
    x3 = states[2]

    # Unwrap states
    u1 = inputs[0]
    u2 = inputs[1]
    u3 = inputs[2]

    # Parameters
    c = np.array([10., 28., 8/3])

    # ODE
    dx1 = c[0]*(x2-x1) + u1
    dx2 = c[1] *x1 - x2 - x1*x3 + u2
    dx3 = x1*x2 - c[2]*x3 + u3


    model.set_dynamical_equations([dx1, dx2, dx3])

    # Initial conditions
    x0 = [0, 0, 0]

    # Create model and run simulation
    dt = 0.001
    model.setup(dt=dt)
    model.set_initial_conditions(x0=x0)
    return model

my attempt to initialize the particle filter ends in the following error:

18 observer= PF(model, roughening=True, plot_backend='bokeh')
---> 19 observer.setup()
20

C:\ProgramData\Anaconda3\lib\site-packages\hilo_mpc\modules\estimator\pf.py in setup(self, **kwargs)
290
291 self._propagate_particles(n_s)
--> 292 self._evaluate_likelihood(n_s)
293
294 n_x = self._model.n_x

C:\ProgramData\Anaconda3\lib\site-packages\hilo_mpc\modules\estimator\pf.py in _evaluate_likelihood(self, n_samples)
156 R = ca.SX.sym('R', (n_y, n_y))
157
--> 158 q = self._normpdf(Y, y, ca.sqrt(R)) # Since R is usually a diagonal matrix, we can use ca.sqrt() here. We need
159 # to change this, if we have non-diagonal matrices, i.e. covariance entries.
160 q /= ca.sum2(q)

C:\ProgramData\Anaconda3\lib\site-packages\casadi\casadi.py in call(self, *args, **kwargs)
23379 if len(args)>0:
23380 # Ordered inputs -> return tuple

23381 ret = self.call(args)
23382 if len(ret)==0:
23383 return None

C:\ProgramData\Anaconda3\lib\site-packages\casadi\casadi.py in call(self, *args)
20037
20038 """

20039 return _casadi.Function_call(self, *args)
20040
20041

RuntimeError: Error in Function::call for 'normpdf' [SXFunction] at .../casadi/core/function.cpp:339:
.../casadi/core/function_internal.hpp:1644: Input 0 (i0) has mismatching shape. Got 3-by-15. Allowed dimensions, in general, are:

• The input dimension N-by-M (here 1-by-1)
• A scalar, i.e. 1-by-1
• M-by-N if N=1 or M=1 (i.e. a transposed vector)
• N-by-M1 if K*M1=M for some K (argument repeated horizontally)
• N-by-P*M, indicating evaluation with multiple arguments (P must be a multiple of 1 for consistency with previous inputs)

Some problem with initialization of probability distribution. I don't understand casadi, so I can't fix it myself.
The code works perfectly with EKF, UKF and MHE filters.
Thanks for help in advance.

Add basic functionalities to LMPC

These are

1. Time-varying parameters (at the moment they do not seem to work)
2. Plotting functionalities ( the solution object needs to be added to the class)
3. Add box constraints

Implement simple Stochastic MPC

Add simple stochastic MPC as in this paper.

Main idea

First of all we should create a stochastic formulation without considering the GP as in the paper.

1. Create new class called something like GaussianRandomVariable to define the stochastic variable. This should contain mean and standard deviation
2. Allow to add (later to multiply) this variable with the model.
3. The SMPC class should detect the different parts (random variables and non) to be able to do the necessary operations.

Note

At the moment one can sum the GP to the RHS of the model using the '+' operator. This adds the mean of the GP to the model. Think about using another operator to add the mean, while keeping the '+' operator to define a stochastic model.

GP Discrete Time State Space Models

@brunomorampc I am trying to use a GP Model in an MPC setting where I am trying to write the state-space equations as

x1(t+1) = gp1(x(t), u(t))
x2(t+1) = gp2(x(t),u(t)),

where x = [x1;x2], u = [u1;u2]

I have currently written code as follows. However I am still stuck at writing the correct syntax for writing gp into the equation. I would really appreciate some help with this.

Thank you in advance!

Create GP Prediction Model

gp_model = Model(plot_backend='bokeh', name='GP_Model', discrete=True)

Set states and inputs

s = gp_model.set_dynamical_states(['x1','x2'])
u = gp_model.set_inputs(['u1','u2'])
gp_model.set_measurement_equations([s[0], s[1]])

x1 = s[0]
x2 = s[1]

u1 = u[0]
u2 = u[1]

x_next = gp1.predict([x1,x2,u1,u2])
y_next = gp2.predict([x1,x2,u1,u2])

gp_model.set_dynamical_equations(['x1', 'x2' ])

Set up the model

gp_model.setup(dt=0.01)

Set initial conditions

x0_GP = [1, 1]
gp_model.set_initial_conditions(x0=x0_GP)

Failed to install the package

I failed to install this package on my computer.

HSL linear solvers

Hi,

Can I use HSL ma27/ ma57 instead of MUMPS for Ipopt using your software ?

Thx

ValueError: Using HILO

Discussed in #26

^{Originally posted by wissamkafa July 31, 2023}
Hi,
First of all, thank you for this wonderful library and the efforts you put into it.
I am interested in using NN to improve traditional controller performance.
I started to use HILO.
I have all packages installed, and following the code written in the documentation:
https://hilo-mpc.github.io/hilo-mpc/examples/nmpc_hybrid_bio.html
I tried to test the first code.

• First in (from hilo_mpc.library.models import ecoli_D1210_conti_plant) the "ecoli_D1210_conti_plant" is not reocgnizable instead I used "ecoli_D1210_conti" --- are they the same? or there is something missing here?
• Second: the code runs without problems until the line "control_loop.run(100, p=[100, 4])"
  here we get an error:
  /home/user/anaconda3/envs/HILO/lib/python3.9/site-packages/hilo_mpc/modules/dynamic_model/dynamic_model.py", line 3965, in simulate
  raise ValueError(f"Dimension mismatch. Supplied dimension for the parameters 'p' is "
  ValueError: Dimension mismatch. Supplied dimension for the parameters 'p' is 2x1, but required dimension is 5x1.

it looks like a dimensions mismatch, and to test this I passed [100, 4, 0, 0, 0] instead [100, 4]
but another error raises:
File "/home/user/anaconda3/envs/HILO/lib/python3.9/site-packages/hilo_mpc/modules/controller/mpc.py", line 763, in optimize
raise ValueError(
ValueError: The model has 2 constant parameter(s): ['Sf', 'If']. You must pass me the value of these before running the optimization to the 'cp' parameter.
I also tried to assign values for Sf and If before calling control_loop.run() but without any improvement.

Any idea why this is happening, am I missing something or what?

Version of numpy

Thank you for a great project! Is it really necessary to border the numpy version numpy<=1.19.5 in setup.cfg?

Bug when using multiple runs in the MPC and one of the run fails

Seems like if one of the runs does not work it returns the failing run instead of the best-working one

parameter varying stage and terminal costs

Hi,
I am wondering if it is possible to implement a parameter varying MPC using hilo?
I mean, suppose that I have a learning model which adjusts the parameters of the MPC at each run. particularly:

1. parameters of the dynamical model, i.e., the equations
2. stage cost parameters
3. terminal cost parameters

That would be perfect if there is an example, you share with me.

Change append to concat in pandas dataframes to be compatible with Pandas >=2.0.0

See issue https://stackoverflow.com/questions/75956209/dataframe-object-has-no-attribute-append

Fix mpc.solution.plot() when model has constant/time-varying parameters

We get a dimension error for this case because the parameters are not added to the solution object at every time steps.

Bokeh plot not working in Jupyter

when running

model.solution.plot()

with the bokeh backend we get an error in Jupyter