Functions: Generic Functions

Functions: Generic Functions

This page contains the documentation of functions that are generic to SimModel, StateEstimator and PredictiveController types.

Set Constraint

# ModelPredictiveControl.setconstraint! — Function

setconstraint!(estim::MovingHorizonEstimator; <keyword arguments>) -> estim

Set the constraint parameters of the MovingHorizonEstimator estim.

It supports both soft and hard constraints on the estimated state , process noise and sensor noise :

and also . All the constraint parameters are vector. Use ±Inf values when there is no bound. The constraint softness parameters , also called equal concern for relaxation, are non-negative values that specify the softness of the associated bound. Use 0.0 values for hard constraints (default for all of them). Notice that constraining the estimated sensor noises is equivalent to bounding the innovation term, since . See Extended Help for details on model augmentation and time-varying constraints.

Arguments

The default constraints are mentioned here for clarity but omitting a keyword argument will not re-assign to its default value (defaults are set at construction only). The same applies for PredictiveController.

estim::MovingHorizonEstimator : moving horizon estimator to set constraints
x̂min=fill(-Inf,nx̂) / x̂max=fill(+Inf,nx̂) : estimated state bound
ŵmin=fill(-Inf,nx̂) / ŵmax=fill(+Inf,nx̂) : estimated process noise bound
v̂min=fill(-Inf,nym) / v̂max=fill(+Inf,nym) : estimated sensor noise bound
c_x̂min=fill(0.0,nx̂) / c_x̂max=fill(0.0,nx̂) : x̂min / x̂max softness weight
c_ŵmin=fill(0.0,nx̂) / c_ŵmax=fill(0.0,nx̂) : ŵmin / ŵmax softness weight
c_v̂min=fill(0.0,nym) / c_v̂max=fill(0.0,nym) : v̂min / v̂max softness weight
all the keyword arguments above but with a first capital letter, e.g. X̂max or C_ŵmax: for time-varying constraints (see Extended Help)

Examples

julia> estim = MovingHorizonEstimator(LinModel(ss(0.5,1,1,0,1)), He=3);

julia> estim = setconstraint!(estim, x̂min=[-50, -50], x̂max=[50, 50])
MovingHorizonEstimator estimator with a sample time Ts = 1.0 s, OSQP optimizer, LinModel and:
 3 estimation steps He
 1 manipulated inputs u (0 integrating states)
 2 estimated states x̂
 1 measured outputs ym (1 integrating states)
 0 unmeasured outputs yu
 0 measured disturbances d

Extended Help

Note that the state and process noise constraints are applied on the augmented model, detailed in SteadyKalmanFilter Extended Help. For variable constraints, the bounds can be modified after calling updatestate!, that is, at runtime, except for ±Inf bounds. Time-varying constraints over the estimation horizon are also possible, mathematically defined as:

For this, use the same keyword arguments as above but with a first capital letter:

X̂min / X̂max / C_x̂min / C_x̂max : constraints (nx̂*(He+1),).
Ŵmin / Ŵmax / C_ŵmin / C_ŵmax : constraints (nx̂*He,).
V̂min / V̂max / C_v̂min / C_v̂max : constraints (nym*He,).

setconstraint!(mpc::PredictiveController; <keyword arguments>) -> mpc

Set the constraint parameters of the PredictiveController mpc.

The predictive controllers support both soft and hard constraints, defined by:

and also . The last line is the terminal constraints applied on the states at the end of the horizon (see Extended Help). See MovingHorizonEstimator constraints for details on bounds and softness parameters . The output and terminal constraints are all soft by default. See Extended Help for time-varying constraints.

Arguments

mpc::PredictiveController : predictive controller to set constraints
umin=fill(-Inf,nu) / umax=fill(+Inf,nu) : manipulated input bound
Δumin=fill(-Inf,nu) / Δumax=fill(+Inf,nu) : manipulated input increment bound
ymin=fill(-Inf,ny) / ymax=fill(+Inf,ny) : predicted output bound
x̂min=fill(-Inf,nx̂) / x̂max=fill(+Inf,nx̂) : terminal constraint bound
c_umin=fill(0.0,nu) / c_umax=fill(0.0,nu) : umin / umax softness weight
c_Δumin=fill(0.0,nu) / c_Δumax=fill(0.0,nu) : Δumin / Δumax softness weight
c_ymin=fill(1.0,ny) / c_ymax=fill(1.0,ny) : ymin / ymax softness weight
c_x̂min=fill(1.0,nx̂) / c_x̂max=fill(1.0,nx̂) : x̂min / x̂max softness weight
all the keyword arguments above but with a first capital letter, except for the terminal constraints, e.g. Ymax or C_Δumin: for time-varying constraints (see Extended Help)

Examples

julia> mpc = LinMPC(setop!(LinModel(tf(3, [30, 1]), 4), uop=[50], yop=[25]));

julia> mpc = setconstraint!(mpc, umin=[0], umax=[100], Δumin=[-10], Δumax=[+10])
LinMPC controller with a sample time Ts = 4.0 s, OSQP optimizer, SteadyKalmanFilter estimator and:
 10 prediction steps Hp
  2 control steps Hc
  1 manipulated inputs u (0 integrating states)
  2 estimated states x̂
  1 measured outputs ym (1 integrating states)
  0 unmeasured outputs yu
  0 measured disturbances d

Extended Help

Terminal constraints provide closed-loop stability guarantees on the nominal plant model. They can render an unfeasible problem however. In practice, a sufficiently large prediction horizon without terminal constraints is typically enough for stability. Note that terminal constraints are applied on the augmented state vector (see SteadyKalmanFilter for details on augmentation).

For variable constraints, the bounds can be modified after calling moveinput!, that is, at runtime, but not the softness parameters . It is not possible to modify ±Inf bounds at runtime.

To keep a variable unconstrained while maintaining the ability to add a constraint later at runtime, set the bound to an absolute value sufficiently large when you create the controller (but different than ±Inf).

It is also possible to specify time-varying constraints over and horizons. In such a case, they are defined by:

For this, use the same keyword arguments as above but with a first capital letter:

Umin / Umax / C_umin / C_umax : constraints (nu*Hp,).
ΔUmin / ΔUmax / C_Δumin / C_Δumax : constraints (nu*Hc,).
Ymin / Ymax / C_ymin / C_ymax : constraints (ny*Hp,).

Evaluate Output y

# ModelPredictiveControl.evaloutput — Function

evaloutput(model::SimModel, d=[]) -> y

Evaluate SimModel outputs y from model.x states and measured disturbances d.

Calling a SimModel object calls this evaloutput method.

Examples

julia> model = setop!(LinModel(tf(2, [10, 1]), 5.0), yop=[20]);

julia> y = evaloutput(model)
1-element Vector{Float64}:
 20.0

evaloutput(estim::StateEstimator, d=[]) -> ŷ

Evaluate StateEstimator outputs ŷ from estim.x̂ states and disturbances d.

Calling a StateEstimator object calls this evaloutput method.

Examples

julia> kf = SteadyKalmanFilter(setop!(LinModel(tf(2, [10, 1]), 5), yop=[20]));

julia> ŷ = evaloutput(kf)
1-element Vector{Float64}:
 20.0

Update State x

# ModelPredictiveControl.updatestate! — Function

updatestate!(model::SimModel, u, d=[]) -> x

Update model.x states with current inputs u and measured disturbances d.

Examples

julia> model = LinModel(ss(1.0, 1.0, 1.0, 0, 1.0));

julia> x = updatestate!(model, [1])
1-element Vector{Float64}:
 1.0

updatestate!(estim::StateEstimator, u, ym, d=[]) -> x̂

Update estim.x̂ estimate with current inputs u, measured outputs ym and dist. d.

The method removes the operating points with remove_op! and call update_estimate!.

Examples

julia> kf = SteadyKalmanFilter(LinModel(ss(0.1, 0.5, 1, 0, 4.0)));

julia> x̂ = updatestate!(kf, [1], [0]) # x̂[2] is the integrator state (nint_ym argument)
2-element Vector{Float64}:
 0.5
 0.0

updatestate!(mpc::PredictiveController, u, ym, d=[]) -> x̂

Call updatestate! on mpc.estim StateEstimator.

Init State x

# ModelPredictiveControl.initstate! — Function

initstate!(model::SimModel, u, d=[]) -> x

Init model.x with manipulated inputs u and measured disturbances d steady-state.

It calls steadystate!(model, u, d):

If model is a LinModel, the method computes the steady-state of current inputs u and measured disturbances d.
Else, model.x is left unchanged. Use setstate! to manually modify it.

Examples

julia> model = LinModel(tf(6, [10, 1]), 2.0);

julia> u = [1]; x = initstate!(model, u); y = round.(evaloutput(model), digits=3)
1-element Vector{Float64}:
 6.0

julia> x ≈ updatestate!(model, u)
true

initstate!(estim::StateEstimator, u, ym, d=[]) -> x̂

Init estim.x̂ states from current inputs u, measured outputs ym and disturbances d.

The method tries to find a good stead-state for the initial esitmate . It removes the operating points with remove_op! and call init_estimate!:

If estim.model is a LinModel, it finds the steady-state of the augmented model using u and d arguments, and uses the ym argument to enforce that . For control applications, this solution produces a bumpless manual to automatic transfer. See init_estimate! for details.
Else, estim.x̂ is left unchanged. Use setstate! to manually modify it.

If applicable, it also sets the error covariance estim.P̂ to estim.P̂0.

Examples

julia> estim = SteadyKalmanFilter(LinModel(tf(3, [10, 1]), 0.5), nint_ym=[2]);

julia> u = [1]; y = [3 - 0.1]; x̂ = round.(initstate!(estim, u, y), digits=3)
3-element Vector{Float64}:
  5.0
  0.0
 -0.1

julia> x̂ ≈ updatestate!(estim, u, y)
true

julia> evaloutput(estim) ≈ y
true

initstate!(mpc::PredictiveController, u, ym, d=[]) -> x̂

Init the states of mpc.estim StateEstimator and warm start mpc.ΔŨ at zero.

Set State x

# ModelPredictiveControl.setstate! — Function

setstate!(model::SimModel, x)

Set model.x states to values specified by x.

setstate!(estim::StateEstimator, x̂)

Set estim.x̂ states to values specified by x̂.

setstate!(mpc::PredictiveController, x̂)

Set the estimate at mpc.estim.x̂.

Quick Simulation

Simulate

# ModelPredictiveControl.sim! — Function

sim!(plant::SimModel, N::Int, u=plant.uop.+1, d=plant.dop; x0=zeros(plant.nx)) -> res

Open-loop simulation of plant for N time steps, default to unit bump test on all inputs.

The manipulated inputs and measured disturbances are held constant at u and d values, respectively. The plant initial state is specified by x0 keyword arguments. The function returns SimResult instances that can be visualized by calling plot from Plots.jl on them (see Examples below). Note that the method mutates plant internal states.

Examples

julia> plant = NonLinModel((x,u,d)->0.1x+u+d, (x,_)->2x, 10.0, 1, 1, 1, 1);

julia> res = sim!(plant, 15, [0], [0], x0=[1])
Simulation results of NonLinModel with 15 time steps.

julia> using Plots; plot(res, plotu=false, plotd=false, plotx=true)

sim!(
    estim::StateEstimator,
    N::Int,
    u = estim.model.uop .+ 1,
    d = estim.model.dop;
    <keyword arguments>
) -> res

Closed-loop simulation of estim estimator for N steps, default to input bumps.

See Arguments for the available options. The noises are provided as standard deviations σ vectors. The simulated sensor and process noises of plant are specified by y_noise and x_noise arguments, respectively.

Arguments

estim::StateEstimator : state estimator to simulate
N::Int : simulation length in time steps
u = estim.model.uop .+ 1 : manipulated input value
d = estim.model.dop : plant measured disturbance value
plant::SimModel = estim.model : simulated plant model
u_step = zeros(plant.nu) : step load disturbance on plant inputs
u_noise = zeros(plant.nu) : gaussian load disturbance on plant inputs
y_step = zeros(plant.ny) : step disturbance on plant outputs
y_noise = zeros(plant.ny) : additive gaussian noise on plant outputs
d_step = zeros(plant.nd) : step on measured disturbances
d_noise = zeros(plant.nd) : additive gaussian noise on measured dist.
x_noise = zeros(plant.nx) : additive gaussian noise on plant states
x0 = plant.x : plant initial state
x̂0 = nothing : initial estimate , initstate! is used if nothing
lastu = plant.uop : last plant input for initialization

Examples

julia> model = LinModel(tf(3, [30, 1]), 0.5);

julia> estim = KalmanFilter(model, σR=[0.5], σQ=[0.25], σQint_ym=[0.01], σP0int_ym=[0.1]);

julia> res = sim!(estim, 50, [0], y_noise=[0.5], x_noise=[0.25], x0=[-10], x̂0=[0, 0])
Simulation results of KalmanFilter with 50 time steps.

julia> using Plots; plot(res, plotŷ=true, plotu=false, plotxwithx̂=true)

sim!(
    mpc::PredictiveController,
    N::Int,
    ry = mpc.estim.model.yop .+ 1,
    d  = mpc.estim.model.dop,
    ru = mpc.estim.model.uop;
    <keyword arguments>
) -> res

Closed-loop simulation of mpc controller for N steps, default to output setpoint bumps.

The output and manipulated input setpoints are held constant at ry and ru, respectively. The keyword arguments are identical to sim!(::StateEstimator, ::Int).

Examples

julia> model = LinModel([tf(3, [30, 1]); tf(2, [5, 1])], 4);

julia> mpc = setconstraint!(LinMPC(model, Mwt=[0, 1], Nwt=[0.01], Hp=30), ymin=[0, -Inf]);

julia> res = sim!(mpc, 25, [0, 0], y_noise=[0.1], y_step=[-10, 0])
Simulation results of LinMPC with 25 time steps.

julia> using Plots; plot(res, plotry=true, plotŷ=true, plotymin=true, plotu=true)

Simulation Results

# ModelPredictiveControl.SimResult — Type

SimResult(
    obj::Union{SimModel, StateEstimator, PredictiveController},
    U_data,
    Y_data,
    D_data  = [];
    X_data  = nothing,
    X̂_data  = nothing,
    Ry_data = nothing,
    Ru_data = nothing
)

Manually construct a SimResult to quickly plot obj simulations.

Except for obj, all the arguments should be matrices of N columns, where N is the number of time steps. SimResult objects allow to quickly plot simulation results. Simply call plot from Plots.jl on them.

Examples

julia> plant = LinModel(tf(1, [1, 1]), 1.0); N = 5; U_data = fill(1.0, 1, N);

julia> Y_data = reduce(hcat, (updatestate!(plant, U_data[:, i]); plant()) for i=1:N)
1×5 Matrix{Float64}:
 0.632121  0.864665  0.950213  0.981684  0.993262

julia> res = SimResult(plant, U_data, Y_data)
Simulation results of LinModel with 5 time steps.

julia> using Plots; plot(res)

Get Additional Information

# ModelPredictiveControl.getinfo — Function

getinfo(estim::MovingHorizonEstimator) -> info

Get additional info on estim MovingHorizonEstimator optimum for troubleshooting.

The function should be called after calling updatestate!. It returns the dictionary info with the following fields:

:Ŵ : optimal estimated process noise over ,
:x̂arr: optimal estimated state at arrival,
:ϵ : optimal slack variable,
:J : objective value optimum,
:X̂ : optimal estimated states over ,
:x̂ : optimal estimated state for the next time step,
:V̂ : optimal estimated sensor noise over ,
:P̄ : estimation error covariance at arrival,
:x̄ : optimal estimation error at arrival,
:Ŷ : optimal estimated outputs over ,
:Ŷm : optimal estimated measured outputs over ,
:Ym : measured outputs over ,
:U : manipulated inputs over ,
:D : measured disturbances over ,
:sol : solution summary of the optimizer for printing

Examples

julia> estim = MovingHorizonEstimator(LinModel(ss(1.0, 1.0, 1.0, 0, 1)), He=1, nint_ym=0);

julia> updatestate!(estim, [0], [1]);

julia> round.(getinfo(estim)[:Ŷ], digits=3)
1-element Vector{Float64}:
 0.5

getinfo(mpc::PredictiveController) -> info

Get additional info about mpc PredictiveController optimum for troubleshooting.

The function should be called after calling moveinput!. It returns the dictionary info with the following fields:

:ΔU : optimal manipulated input increments over ,
:ϵ : optimal slack variable,
:J : objective value optimum,
:U : optimal manipulated inputs over ,
:u : current optimal manipulated input,
:d : current measured disturbance,
:D̂ : predicted measured disturbances over ,
:ŷ : current estimated output,
:Ŷ : optimal predicted outputs over ,
:x̂end: optimal terminal states,
:Ŷs : predicted stochastic output over of InternalModel,
:R̂y : predicted output setpoint over ,
:R̂u : predicted manipulated input setpoint over ,

For LinMPC and NonLinMPC, the field :sol also contains the optimizer solution summary that can be printed. Lastly, the optimal economic cost :JE is also available for NonLinMPC.

Examples

julia> mpc = LinMPC(LinModel(tf(5, [2, 1]), 3), Nwt=[0], Hp=1, Hc=1);

julia> u = moveinput!(mpc, [10]);

julia> round.(getinfo(mpc)[:Ŷ], digits=3)
1-element Vector{Float64}:
 10.0