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 boundary conditions. The OCP thus requires a prediction of these boundary conditions 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 boundary conditions 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.boundaryconditions 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 emulator, control and boundaryconditions objects are defined, an mpc object can be instantiated:

mpc = mpcpy.MPC(emulator,control,boundaryconditions,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 boundary conditions 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¶

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

# instantiate the emulator
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.)
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)

Boundary conditions¶

The mpcpy.boundaryconditions class is a helper class for manging the required boundary conditions. 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.boundaryconditions 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 = np.random.random(len(time)) - 0.5
boundaryconditions = mpcpy.Boundaryconditions({'time':time,'d':d})

bcs = boundaryconditions(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 boundary conditions 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 contol inputs
        optimizer = pyomo.SolverFactory('cplex')
        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.boundaryconditions object which represent the actual conditions. From these values predictions may be derived.

Example¶

prediction = mpcpy.Prediction(boundaryconditions)
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,boundaryconditions,emulationtime=1000,resulttimestep=10)

# run the mpc
res = mpc()
print(res)