parameter varying stage and terminal costs

[aabsz]

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.

[brunomorampc]

Hi @aabsz , yes there is a way. It a bit depends on what you want to do. If you already have the time-varying profiles of the parameters you can define them as time-varying parameters (see option 1). If you want to change them online then you can just pass the new value online (see option 2). Note that the two approaches are conceptually different but the may lead to similar (not the same!) results. This section of the documentation explains how time-varying parameters are defined.

Option 1: you want to pass the time-varying parameter profiles

Time-varying model parameters

The trick is defining the parameters of the model, then set these parameters as "time-varying" in controller. As an example, have a look as this pendulum subject to a sinusoidal disturbance:

from hilo_mpc import Model, NMPC
import casadi as ca
import numpy as np
model = Model(plot_backend='bokeh')
# Constants
M = 5.
m = 1.
l = 1.
h = .5
g = 9.81

# States and algebraic variables
x = model.set_dynamical_states(['x', 'v', 'theta', 'omega'])
model.set_measurements(['yx', 'yv', 'ytheta', 'tomega'])
model.set_measurement_equations([x[0], x[1], x[2], x[3]])
d = model.set_parameters(['d_tvp', 'd_const'])  # disturbance
v = x[1]
theta = x[2]
omega = x[3]
# Inputs
F = model.set_inputs('F')

# ODEs

dx = v
dv = 1. / (M + m - m * ca.cos(theta)) * (m * g * ca.sin(theta) - m * l * ca.sin(theta) * omega ** 2 + F)
dtheta = omega + d[0]
domega = 1. / l * (dv * ca.cos(theta) + g * ca.sin(theta)) + d[1]

# Algebraic equations
# dy = h + l * ca.cos(theta) - y

model.set_equations(ode=[dx, dv, dtheta, domega])

# Initial conditions
x0 = [2.5, 0., 0.1, 0.]
u0 = 0.

# Create model and run simulation
dt = .1
model.setup(dt=dt)

nmpc = NMPC(model)
nmpc.quad_stage_cost.add_states(names=['v', 'theta'], ref=[0, 0], weights=[10, 5])
nmpc.quad_stage_cost.add_inputs(names='F', weights=0.1)
nmpc.horizon = 25
nmpc.set_box_constraints(x_ub=[5, 10, 10, 10], x_lb=[-5, -10, -10, -10])
nmpc.set_initial_guess(x_guess=x0, u_guess=u0)
nmpc.set_time_varying_parameters(names=['d_tvp'], values={'d_tvp': np.sin(np.arange(0, 50, dt)).tolist()})
nmpc.setup()

n_steps = 50
model.set_initial_conditions(x0=x0)
sol = model.solution

for step in range(n_steps):
    u = nmpc.optimize(x0, cp=0.1)
    model.simulate(u=u, p=[ca.sin(dt * step), 0.1])
    x0 = sol['x:f']

sol.plot()

Time-varying weights of the MPC

This works also for the weights of the controller: the trick is to use the parameters of the model as weights for the MPC. This is sort of a work around, because defining the parameter of the controller in the model is something that we wouldn't like to do, but it works: For example, this point mass system

from hilo_mpc import Model, NMPC

# Define the model
model = Model(plot_backend='bokeh')

# Constants
M = 5.

# States and algebraic variables
xx = model.set_dynamical_states('x', 'vx', 'y', 'vy')
p = model.set_parameters('w_x', 'w_y')
x = xx[0]
vx = xx[1]
y = xx[2]
vy = xx[3]

# Inputs
F = model.set_inputs('Fx', 'Fy')
Fx = F[0]
Fy = F[1]

# ODEs
dd1 = vx
dd2 = Fx / M
dd3 = vy
dd4 = Fy / M

model.set_dynamical_equations([dd1, dd2, dd3, dd4])

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

# time interval
dt = 0.1
model.setup(dt=dt)
model.set_initial_conditions(x0=x0)
sol = model.solution
n_steps = 200
nmpc = NMPC(model, stats=True)
nmpc.stage_cost.cost = p[0] * (x - 2) ** 2 + p[1] * (y - 3) ** 2
nmpc.horizon = 20
nmpc.quad_stage_cost.add_inputs(names=['Fx', 'Fy'], weights=[0.1, 0.1])
nmpc.set_box_constraints(x_ub=[50, 10, 10, 10], x_lb=[-10, -10, -10, -10])
nmpc.set_initial_guess(x_guess=x0, u_guess=u0)
tvp = {'w_x': [0] * int(n_steps / 2) + [10] * int(n_steps / 2),
       'w_y': [10] * int(n_steps / 2) + [0] * int(n_steps / 2)}
nmpc.set_time_varying_parameters(names=['w_x', 'w_y'],
                                 values=tvp)

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

for step in range(n_steps):
    u = nmpc.optimize(x0)
    model.simulate(u=u, p=[tvp['w_x'][step], tvp['w_y'][step]])
    x0 = sol['x:f']

model.solution.plot()

Note that model.simulate expects 2 parameters, since they were defined in the model. They have of course no effect in the dynamic but you still need to pass them. One could pass simply any value, but if you want to see the correct values in the plots, then make sure to pass the value of the parameter at the correct time step.

Option 2: you want to change the parameters/weights online

Then it is sufficient to proceed, as before, defining the parameters of the controller/model, but you do not need to define them as time-varying. You would need only to supply the value of the parameters in the for loop as follows:

"""
In this case the weights are not defined as tvp but only changed online
"""

from hilo_mpc import Model, NMPC

# Define the model
model = Model(plot_backend='bokeh')

# Constants
M = 5.

# States and algebraic variables
xx = model.set_dynamical_states('x', 'vx', 'y', 'vy')
p = model.set_parameters('w_x', 'w_y')
x = xx[0]
vx = xx[1]
y = xx[2]
vy = xx[3]

# Inputs
F = model.set_inputs('Fx', 'Fy')
Fx = F[0]
Fy = F[1]

# ODEs
dd1 = vx
dd2 = Fx / M
dd3 = vy
dd4 = Fy / M

model.set_dynamical_equations([dd1, dd2, dd3, dd4])

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

# time interval
dt = 0.1
model.setup(dt=dt)
model.set_initial_conditions(x0=x0)
sol = model.solution
n_steps = 200
nmpc = NMPC(model, stats=True)
nmpc.stage_cost.cost = p[0] * (x - 2) ** 2 + p[1] * (y - 3) ** 2
nmpc.horizon = 20
nmpc.quad_stage_cost.add_inputs(names=['Fx', 'Fy'], weights=[0.1, 0.1])
nmpc.set_box_constraints(x_ub=[50, 10, 10, 10], x_lb=[-10, -10, -10, -10])
nmpc.set_initial_guess(x_guess=x0, u_guess=u0)
w_x = [0] * int(n_steps / 2) + [10] * int(n_steps / 2)
w_y = [0] * int(n_steps / 2) + [10] * int(n_steps / 2)
nmpc.setup()
for step in range(n_steps):
    u = nmpc.optimize(x0, cp = [w_x[step], w_y[step]])
    model.simulate(u=u, p=[w_x[step], w_y[step]])
    x0 = sol['x:f']

model.solution.plot()

Please let me know if something is not clear or you have further questions

[aabsz]

Dear Bruno,

Thank you so much for the comprehensive reply.

I have a question regarding "option 1: Time-varying model parameters."

In your example, you have defined d[0] as a variable parameter of the dynamical model. I noticed that before setting up the NMPC model, you defined the variable parameter as follows:

nmpc.set_time_varying_parameters(names=['d_tvp'], values={'d_tvp': np.sin(np.arange(0, 50, dt)).tolist()})

Then, when creating the control loop, you specified its exact value like this:

model.simulate(u=u, p=[np.sin(dt * step), 0.1]) (I assume "p" is incorrect, and it should be "d." Am I correct?)

In this context, you assumed that you know the values of this variable parameter beforehand. Now, my question is, what if I don't know its value in advance? How can I pass it in the simulation? Suppose that its value is calculated somewhere in my code, but I don't know how to define it when setting up the NMPC. In other words, I'm uncertain about how to define it in the following line:

nmpc.set_time_varying_parameters(names=['d_tvp'], values={'d_tvp': np.sin(np.arange(0, 50, dt)).tolist()})

[aabsz]

Hi, I figured it out!