Quickstart¶
MPC programming structure¶
In Model Predictive Control (MPC) a dynamic system is controlled with some control inputs. These inputs are computed with the aid of a model of the dynamic system which predicts the system behavior over a certain time horizon.
An MPC simulation is here defined as a dynamic simulation of such a system. In such a simulation the actual system is replaced by a system emulator, a dynamic simulation which should model the system behavior a close as possible.
The control algorithm contains a more simplified model of the dynamic system and often uses an optimization technique to compute the control signals which minimize some objective function over a time horizon. This is referred to as an Optimal Control Problem (OCP).
The system under consideration is most likely subjected to some uncontrolled disturbances. The OCP thus requires a prediction of these disturbances over the prediction horizon to be able to compute the system behavior. Furthermore, an estimation of the current system state is needed, probably requiring some (virtual) measurements.
An MPC simulation can thus be broken down into the following steps:
- Estimate the state of the system.
- Generate a prediction for the disturbances over the control horizon.
- Compute control signals by solving an OCP.
- Supply the control signals to a system emulator and simulate the system behavior.
- Repeat.
In mpcpy the mpcpy.mcp object contains the main functions for
orchestrating the MPC. It depends on an mpcpy.emulator object, an
mpcpy.disturbances object and an mpcpy.control object.
The emulation time is split up in parts defined by the control receding time. For each part, the above steps are executed and the emulator simulates the system over the receding time.
The mpcpy.control object handles the generation of control signals.
Therfore it requires an mpcpy.stateestimation object and an
mpcpy.prediction object.
Example¶
When the mpcpy.emulator, mpcpy.control and
mpcpy.disturbances objects are defined, an mpcpy.mpc object can
be instantiated:
mpc = mpcpy.MPC(emulator, control, disturbances, emulationtime=1000, resulttimestep=10)
Emulator¶
The mpcpy.emulator class serves as a base class for defining emulator
objects. By default, the emulator object returns the control solution as its
result. So when no detailed emulator model is available, the
mpcpy.emulator class can be used to split an OCP in different parts and
solve it sequentially. Care must be taken that the OCPs overlap sufficiently, so
the receding time must be significantly shorter than the horizon.
The mpcpy.emulator class requires a list of keys as input arguments,
this list specifies which disturbances and control signals should be
passed to the simulation, all control signals are passed to the simulation.
When subclassing the mpcpy.emulator class the main method to be
redefined is the simulate method. This method receives a starttime,
stoptime and input dictionary as arguments and should return a results
dictionary with a time array and values for the results between the start- and
stoptime.
Internally a res dictionary holds the results of the simulation up to
the current point in time. This res dictionary can also be used to
determine the state required for the OCP.
The system can then be simulated by calling an mpcpy.emulator object
with an array of input times (the time steps at which results are required) and
a dictionary of inputs. During an MPC simulation this is handled by the
mpcpy.mpc object.
Example¶
# Define an emulator class
class Emulator(mpcpy.Emulator):
"""
A custom system emulator
"""
def simulate(self, starttime, stoptime, input):
dt = 1
time = np.arange(starttime, stoptime+dt, dt, dtype=np.float)
# initialize
x = np.ones_like(time)*self.res['x'][-1]
# interpolate inputs
u = np.interp(time, input['time'], input['u'])
d = np.interp(time, input['time'], input['d'])
# perform simulation
for i, t in enumerate(time[:-1]):
# dx/dt = A*x + d + u
x[i+1] = x[i] + (self.parameters['A']*x[i] + d[i] + u[i])*dt
# create and return a results dict
res = {
'time': time,
'x': x,
'd': d,
'u': u,
}
return res
emulator = Emulator(['u', 'd'], parameters={'A': -0.2}, initial_conditions={'x': 0})
# test the emulator with some random data
time = np.arange(0., 1001., 10.)
np.random.seed(0)
d = np.random.random(len(time)) - 0.5
u = 1.0*np.ones(len(time))
emulator.initialize()
res = emulator(time, {'time': time, 'd': d, 'u': u})
print(res)
Disturbances¶
The mpcpy.disturbances class is a helper class for manging the
disturbances.
Signals are defined based on a dictionary which must include a time key
and specifies the values of all signals at these time steps.
An mpcpy.disturbances object can then interpolate the values to
the required time points in the mpcpy.emulator and
mpcpy.prediction objects.
Example¶
time = np.arange(0., 1001., 10.)
d = 0.5*np.sin(2*np.pi*time/1000)
disturbances = mpcpy.Disturbances({'time': time, 'd': d})
bcs = disturbances(np.array([0, 20, 40, 60, 100]))
print(bcs)
Control¶
The emulator requires control inputs which are generated by an object derived
from the mpcpy.Control base class.
The methods that should be redefined in a child class are the
formulation and solution methods.
The formulation method is called once by the mpcpy.mpc object
before the start of the simulation. This method can thus be used to formulate an
abstract optimal control problem and assign values that will not change through
the course of the simulation.
The :code`solution` method should generate the control signals to be used by the
emulator. When the control signals should be generated, it is called with two
parameters representing the current state of the system and the predictions of
the disturbances over the control horizon. These are generated by the
mpcpy.stateestimation and mpcpy.prediction objects supplied to
the control object during initialization respectively.
The :code`solution` method must return a dictionary with a time vector and the
control signal values over the entire control horizon.
Example¶
class SetpointControl(mpcpy.Control):
"""
A control to keep the state as close to a set point as possible
"""
def formulation(self):
# create a pyomo model
model = pyomo.AbstractModel()
model.i = pyomo.Set()
model.ip = pyomo.Set()
model.time = pyomo.Param(model.ip)
model.d = pyomo.Param(model.ip, initialize=0.)
model.x = pyomo.Var(model.ip, domain=pyomo.Reals, initialize=0.)
model.u = pyomo.Var(model.ip, domain=pyomo.NonNegativeReals, bounds=(0., 1.), initialize=0.)
model.x0 = pyomo.Param(initialize=0.)
model.initialcondition = pyomo.Constraint(
rule=lambda model: model.x[0] == model.x0
)
model.constraint = pyomo.Constraint(
model.i,
rule=lambda model, i: (model.x[i+1]-model.x[i])/(model.time[i+1]-model.time[i]) ==
self.parameters['A']*model.x[i] + model.d[i] + model.u[i]
)
model.objective = pyomo.Objective(
rule=lambda model: sum((model.x[i]-self.parameters['set'])**2 for i in model.i)
)
# store the model inside the object
self.model = model
def solution(self, sta, pre):
# create data and instantiate the pyomo model
ip = np.arange(len(pre['time']))
data = {
None: {
'i': {None: ip[:-1]},
'ip': {None: ip},
'time': {(i,): v for i, v in enumerate(pre['time'])},
'x0': {None: sta['x']},
'd': {(i,): pre['d'][i] for i in ip},
}
}
instance = self.model.create_instance(data)
# solve and return the control inputs
optimizer = pyomo.SolverFactory('ipopt')
results = optimizer.solve(instance)
sol = {
'time': np.array([pyomo.value(instance.time[i]) for i in instance.ip]),
'x': np.array([pyomo.value(instance.x[i]) for i in instance.ip]),
'u': np.array([pyomo.value(instance.u[i]) for i in instance.ip]),
'd': np.array([pyomo.value(instance.d[i]) for i in instance.ip]),
}
return sol
When the stateestimation and prediction objects are defined, a control object can be instantiated:
control = SetpointControl(stateestimation, prediction, parameters={'A': -0.2, 'set': 3.0},
horizon=100., timestep=10., receding=10.)
State estimation¶
The mpcpy.stateestimation class is a base class for defining object that
estimate the current state required for the control based on measurements done
on the emulator. When the states required for the control problem are not
directly measurable they are often computed with a Kalman filter.
The main method to redefine in a child class is the :code`stateestimation`
method. This method takes the time at which the states should be estimated as an
input parameter and must return a dictionary of key-value pairs representing the
state at that time.
Example¶
class StateestimationPerfect(mpcpy.Stateestimation):
"""
Perfect state estimation
"""
def stateestimation(self, time):
return {'x': np.interp(time, self.emulator.res['time'], self.emulator.res['x'])}
stateestimation = StateestimationPerfect(emulator)
sta = stateestimation(0)
print(sta)
Prediction¶
The mpcpy.prediction class is a base class for returning prediction data
required for the OCP.
It requires a mpcpy.disturbances object which represent the actual
conditions. From these values predictions may be derived.
Example¶
prediction = mpcpy.Prediction(disturbances)
pre = prediction(np.array([0, 20, 40, 60, 100]))
print(pre)
MPC¶
The entire MPC simulation can now be run from the above defined objects. Just
call the mpcpy.MPC object to start the simulation. Results will be
returned in a dictionary.
Example¶
control = SetpointControl(stateestimation, prediction, parameters={'A': -0.2, 'set': 3.0},
horizon=100., timestep=10., receding=10.)
mpc = mpcpy.MPC(emulator, control, disturbances, emulationtime=1000, resulttimestep=10)
# run the mpc
res = mpc(verbose=1)
print(res)
The optimal control inputs and resulting state of the quickstart example is presented below.