Functions: Generic Functions
This page contains the documentation of functions that are generic to SimModel
, StateEstimator
and PredictiveController
types.
Set Constraint
ModelPredictiveControl.setconstraint!
— Functionsetconstraint!(estim::MovingHorizonEstimator; <keyword arguments>) -> estim
Set the bound constraint parameters of the MovingHorizonEstimator
estim
.
It supports both soft and hard constraints on the estimated state $\mathbf{x̂}$, process noise $\mathbf{ŵ}$ and sensor noise $\mathbf{v̂}$:
\[\begin{alignat*}{3} \mathbf{x̂_{min} - c_{x̂_{min}}} ϵ ≤&&\ \mathbf{x̂}_k(k-j+p) &≤ \mathbf{x̂_{max} + c_{x̂_{max}}} ϵ &&\qquad j = N_k, N_k - 1, ... , 0 \\ \mathbf{ŵ_{min} - c_{ŵ_{min}}} ϵ ≤&&\ \mathbf{ŵ}(k-j+p) &≤ \mathbf{ŵ_{max} + c_{ŵ_{max}}} ϵ &&\qquad j = N_k, N_k - 1, ... , 1 \\ \mathbf{v̂_{min} - c_{v̂_{min}}} ϵ ≤&&\ \mathbf{v̂}(k-j+1) &≤ \mathbf{v̂_{max} + c_{v̂_{max}}} ϵ &&\qquad j = N_k, N_k - 1, ... , 1 \end{alignat*}\]
and also $ϵ ≥ 0$. All the constraint parameters are vector. Use ±Inf
values when there is no bound. The constraint softness parameters $\mathbf{c}$, 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 $\mathbf{v̂}(k) = \mathbf{y^m}(k) - \mathbf{ŷ^m}(k)$. See Extended Help for details on the constant $p$, on model augmentation and on time-varying constraints.
Arguments
All the keyword arguments have non-Unicode alternatives e.g. xhatmin
or Vhatmax
.
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).
estim::MovingHorizonEstimator
: moving horizon estimator to set constraintsx̂min=fill(-Inf,nx̂)
/x̂max=fill(+Inf,nx̂)
: estimated state bound $\mathbf{x̂_{min/max}}$ŵmin=fill(-Inf,nx̂)
/ŵmax=fill(+Inf,nx̂)
: estimated process noise bound $\mathbf{ŵ_{min/max}}$v̂min=fill(-Inf,nym)
/v̂max=fill(+Inf,nym)
: estimated sensor noise bound $\mathbf{v̂_{min/max}}$c_x̂min=fill(0.0,nx̂)
/c_x̂max=fill(0.0,nx̂)
:x̂min
/x̂max
softness weight $\mathbf{c_{x̂_{min/max}}}$c_ŵmin=fill(0.0,nx̂)
/c_ŵmax=fill(0.0,nx̂)
:ŵmin
/ŵmax
softness weight $\mathbf{c_{ŵ_{min/max}}}$c_v̂min=fill(0.0,nym)
/c_v̂max=fill(0.0,nym)
:v̂min
/v̂max
softness weight $\mathbf{c_{v̂_{min/max}}}$- all the keyword arguments above but with a first capital letter, e.g.
X̂max
orC_ŵ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
0 slack variable ϵ (estimation constraints)
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
Extended Help
The constant $p=0$ if estim.direct==true
(current form), else $p=1$ (prediction form). Note that the state $\mathbf{x̂}$ and process noise $\mathbf{ŵ}$ 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 $H_e$ are also possible, mathematically defined as:
\[\begin{alignat*}{3} \mathbf{X̂_{min} - C_{x̂_{min}}} ϵ ≤&&\ \mathbf{X̂} &≤ \mathbf{X̂_{max} + C_{x̂_{max}}} ϵ \\ \mathbf{Ŵ_{min} - C_{ŵ_{min}}} ϵ ≤&&\ \mathbf{Ŵ} &≤ \mathbf{Ŵ_{max} + C_{ŵ_{max}}} ϵ \\ \mathbf{V̂_{min} - C_{v̂_{min}}} ϵ ≤&&\ \mathbf{V̂} &≤ \mathbf{V̂_{max} + C_{v̂_{max}}} ϵ \end{alignat*}\]
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
: $\mathbf{X̂}$ constraints(nx̂*(He+1),)
.Ŵmin
/Ŵmax
/C_ŵmin
/C_ŵmax
: $\mathbf{Ŵ}$ constraints(nx̂*He,)
.V̂min
/V̂max
/C_v̂min
/C_v̂max
: $\mathbf{V̂}$ constraints(nym*He,)
.
setconstraint!(mpc::PredictiveController; <keyword arguments>) -> mpc
Set the bound constraint parameters of the PredictiveController
mpc
.
The predictive controllers support both soft and hard constraints, defined by:
\[\begin{alignat*}{3} \mathbf{u_{min} - c_{u_{min}}} ϵ ≤&&\ \mathbf{u}(k+j) &≤ \mathbf{u_{max} + c_{u_{max}}} ϵ &&\qquad j = 0, 1 ,..., H_p - 1 \\ \mathbf{Δu_{min} - c_{Δu_{min}}} ϵ ≤&&\ \mathbf{Δu}(k+j) &≤ \mathbf{Δu_{max} + c_{Δu_{max}}} ϵ &&\qquad j = 0, 1 ,..., H_c - 1 \\ \mathbf{y_{min} - c_{y_{min}}} ϵ ≤&&\ \mathbf{ŷ}(k+j) &≤ \mathbf{y_{max} + c_{y_{max}}} ϵ &&\qquad j = 1, 2 ,..., H_p \\ \mathbf{x̂_{min} - c_{x̂_{min}}} ϵ ≤&&\ \mathbf{x̂}_i(k+j) &≤ \mathbf{x̂_{max} + c_{x̂_{max}}} ϵ &&\qquad j = H_p \end{alignat*}\]
and also $ϵ ≥ 0$. 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 $\mathbf{c}$. The output and terminal constraints are all soft by default. See Extended Help for time-varying constraints.
Arguments
The keyword arguments Δumin
, Δumax
, c_Δumin
, c_Δumax
, x̂min
, x̂max
, c_x̂min
, c_x̂max
and their capital letter versions have non-Unicode alternatives e.g. Deltaumin
, xhatmax
and C_Deltaumin
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).
mpc::PredictiveController
: predictive controller to set constraintsumin=fill(-Inf,nu)
/umax=fill(+Inf,nu)
: manipulated input bound $\mathbf{u_{min/max}}$Δumin=fill(-Inf,nu)
/Δumax=fill(+Inf,nu)
: manipulated input increment bound $\mathbf{Δu_{min/max}}$ymin=fill(-Inf,ny)
/ymax=fill(+Inf,ny)
: predicted output bound $\mathbf{y_{min/max}}$x̂min=fill(-Inf,nx̂)
/x̂max=fill(+Inf,nx̂)
: terminal constraint bound $\mathbf{x̂_{min/max}}$c_umin=fill(0.0,nu)
/c_umax=fill(0.0,nu)
:umin
/umax
softness weight $\mathbf{c_{u_{min/max}}}$c_Δumin=fill(0.0,nu)
/c_Δumax=fill(0.0,nu)
:Δumin
/Δumax
softness weight $\mathbf{c_{Δu_{min/max}}}$c_ymin=fill(1.0,ny)
/c_ymax=fill(1.0,ny)
:ymin
/ymax
softness weight $\mathbf{c_{y_{min/max}}}$c_x̂min=fill(1.0,nx̂)
/c_x̂max=fill(1.0,nx̂)
:x̂min
/x̂max
softness weight $\mathbf{c_{x̂_{min/max}}}$- all the keyword arguments above but with a first capital letter, except for the terminal constraints, e.g.
Ymax
orC_Δ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 slack variable ϵ (control constraints)
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
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 $H_p$ without terminal constraints is typically enough for stability. If mpc.estim.direct==true
, the estimator computes the states at $i = k$ (the current time step), otherwise at $i = k - 1$. Note that terminal constraints are applied on the augmented state vector $\mathbf{x̂}$ (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 $\mathbf{c}$. 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 $H_p$ and $H_c$ horizons. In such a case, they are defined by:
\[\begin{alignat*}{3} \mathbf{U_{min} - C_{u_{min}}} ϵ ≤&&\ \mathbf{U} &≤ \mathbf{U_{max} + C_{u_{max}}} ϵ \\ \mathbf{ΔU_{min} - C_{Δu_{min}}} ϵ ≤&&\ \mathbf{ΔU} &≤ \mathbf{ΔU_{max} + C_{Δu_{max}}} ϵ \\ \mathbf{Y_{min} - C_{y_{min}}} ϵ ≤&&\ \mathbf{Ŷ} &≤ \mathbf{Y_{max} + C_{y_{max}}} ϵ \end{alignat*}\]
For this, use the same keyword arguments as above but with a first capital letter:
Umin
/Umax
/C_umin
/C_umax
: $\mathbf{U}$ constraints(nu*Hp,)
.ΔUmin
/ΔUmax
/C_Δumin
/C_Δumax
: $\mathbf{ΔU}$ constraints(nu*Hc,)
.Ymin
/Ymax
/C_ymin
/C_ymax
: $\mathbf{Ŷ}$ constraints(ny*Hp,)
.
Evaluate Output y
ModelPredictiveControl.evaloutput
— Functionevaloutput(model::SimModel, d=[]) -> y
Evaluate SimModel
outputs y
from model.x0
states and measured disturbances d
.
It returns model
output at the current time step $\mathbf{y}(k)$. 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=[]) -> ŷ
Evaluate StateEstimator
outputs ŷ
from estim.x̂0
states and disturbances d
.
It returns estim
output at the current time step $\mathbf{ŷ}(k)$. If estim.direct
is true
, the method preparestate!
should be called beforehand to correct the state estimate.
Calling a StateEstimator
object calls this evaloutput
method.
Examples
julia> kf = SteadyKalmanFilter(setop!(LinModel(tf(2, [10, 1]), 5), yop=[20]), direct=false);
julia> ŷ = evaloutput(kf)
1-element Vector{Float64}:
20.0
Change State x
Prepare State x
ModelPredictiveControl.preparestate!
— Functionpreparestate!(model::SimModel, _ , _ ) -> x
Do nothing for SimModel
and return the current model state $\mathbf{x}(k)$.
preparestate!(estim::StateEstimator, ym, d=estim.model.dop) -> x̂
Prepare estim.x̂0
estimate with meas. outputs ym
and dist. d
for the current time step.
This function should be called at the beginning of each discrete time step. Its behavior depends if estim
is a StateEstimator
in the current/filter (1.) or delayed/predictor (2.) formulation:
- If
estim.direct
istrue
, it removes the operating points withremove_op!
, callscorrect_estimate!
, and returns the corrected state estimate $\mathbf{x̂}_k(k)$. - Else, it does nothing and returns the current best estimate $\mathbf{x̂}_{k-1}(k)$.
Examples
julia> estim2 = SteadyKalmanFilter(LinModel(ss(0.1, 0.5, 1, 0, 4)), nint_ym=0, direct=true);
julia> x̂ = round.(preparestate!(estim2, [1]), digits=3)
1-element Vector{Float64}:
0.01
julia> estim1 = SteadyKalmanFilter(LinModel(ss(0.1, 0.5, 1, 0, 4)), nint_ym=0, direct=false);
julia> x̂ = preparestate!(estim1, [1])
1-element Vector{Float64}:
0.0
preparestate!(mpc::PredictiveController, ym, d=[]) -> x̂
Call preparestate!
on mpc.estim
StateEstimator
.
Update State x
ModelPredictiveControl.updatestate!
— Functionupdatestate!(model::SimModel, u, d=[]) -> xnext
Update model.x0
states with current inputs u
and meas. dist. d
for the next time step.
The method computes and returns the model state for the next time step $\mathbf{x}(k+1)$.
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̂next
Update estim.x̂0
estimate with current inputs u
, measured outputs ym
and dist. d
.
This function should be called at the end of each discrete time step. It removes the operating points with remove_op!
, calls update_estimate!
and returns the state estimate for the next time step $\mathbf{x̂}_k(k+1)$. The method preparestate!
should be called prior to this one to correct the estimate when applicable (if estim.direct == true
). Note that the MovingHorizonEstimator
with the default direct=true
option is not able to estimate $\mathbf{x̂}_k(k+1)$, the returned value is therefore the current corrected state $\mathbf{x̂}_k(k)$.
Examples
julia> kf = SteadyKalmanFilter(LinModel(ss(0.1, 0.5, 1, 0, 4.0))); u = [1]; ym = [0];
julia> preparestate!(kf, ym);
julia> x̂ = updatestate!(kf, u, ym) # x̂[2] is the integrator state (nint_ym argument)
2-element Vector{Float64}:
0.5
0.0
updatestate!(mpc::PredictiveController, u, ym, d=[]) -> x̂next
Call updatestate!
on mpc.estim
StateEstimator
.
Init State x
ModelPredictiveControl.initstate!
— Functioninitstate!(model::SimModel, u, d=[]) -> x
Init model.x0
with manipulated inputs u
and meas. dist. d
steady-state.
The method tries to initialize the model state $\mathbf{x}$ at steady-state. It removes the operating points on u
and d
and calls steadystate!
:
- If
model
is aLinModel
, the method computes the steady-state of current inputsu
and measured disturbancesd
. - Else,
model.x0
is left unchanged. Usesetstate!
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̂0
states from current inputs u
, measured outputs ym
and disturbances d
.
The method tries to find a good steady-state for the initial estimate $\mathbf{x̂}$. It stores u - estim.model.uop
at estim.lastu0
and removes the operating points with remove_op!
, and call init_estimate!
:
- If
estim.model
is aLinModel
, it finds the steady-state of the augmented model usingu
andd
arguments, and uses theym
argument to enforce that $\mathbf{ŷ^m}(0) = \mathbf{y^m}(0)$. For control applications, this solution produces a bumpless manual to automatic transfer. Seeinit_estimate!
for details. - Else,
estim.x̂0
is left unchanged. Usesetstate!
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], direct=false);
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.ΔŨ
at zero.
Set State x
ModelPredictiveControl.setstate!
— Functionsetstate!(model::SimModel, x) -> model
Set model.x0
to x - model.xop
from the argument x
.
setstate!(estim::StateEstimator, x̂) -> estim
Set estim.x̂0
to x̂ - estim.x̂op
from the argument x̂
.
setstate!(mpc::PredictiveController, x̂) -> mpc
Set mpc.estim.x̂0
to x̂ - estim.x̂op
from the argument x̂
.
Set Model and Weights
ModelPredictiveControl.setmodel!
— Functionsetmodel!(estim::StateEstimator, model=estim.model; <keyword arguments>) -> estim
Set model
and covariance matrices of estim
StateEstimator
.
Allows model adaptation of estimators based on LinModel
at runtime. Modification of NonLinModel
state-space functions is not supported. New covariance matrices can be specified with the keyword arguments (see SteadyKalmanFilter
documentation for the nomenclature). Not supported by Luenberger
and SteadyKalmanFilter
, use the time-varying KalmanFilter
instead. The MovingHorizonEstimator
model is kept constant over the estimation horizon $H_e$. The matrix dimensions and sample time must stay the same. Note that the observability and controllability of the new augmented model is not verified (see Extended Help for more info).
Arguments
Keyword arguments with emphasis
are non-Unicode alternatives.
estim::StateEstimator
: estimator to set model and covariances.model=estim.model
: new plant model (not supported byNonLinModel
).Q̂=nothing
orQhat
: new augmented model $\mathbf{Q̂}$ covariance matrix.R̂=nothing
orRhat
: new augmented model $\mathbf{R̂}$ covariance matrix.
Examples
julia> kf = KalmanFilter(LinModel(ss(0.1, 0.5, 1, 0, 4.0)), σQ=[√4.0], σQint_ym=[√0.25]);
julia> kf.model.A[], kf.Q̂[1, 1], kf.Q̂[2, 2]
(0.1, 4.0, 0.25)
julia> setmodel!(kf, LinModel(ss(0.42, 0.5, 1, 0, 4.0)), Q̂=[1 0;0 0.5]);
julia> kf.model.A[], kf.Q̂[1, 1], kf.Q̂[2, 2]
(0.42, 1.0, 0.5)
Extended Help
Extended Help
Using the default model augmentation computed by the default_nint
method, switching from a non-integrating plant model to an integrating one will produce an augmented model that is not observable. Moving the unmeasured disturbances at the model input (nint_u
parameter) can fix this issue.
Throw an error if setmodel!
is called on Luenberger
observer w/o the default values.
setmodel!(mpc::PredictiveController, model=mpc.estim.model; <keyword arguments>) -> mpc
Set model
and objective function weights of mpc
PredictiveController
.
Allows model adaptation of controllers based on LinModel
at runtime. Modification of NonLinModel
state-space functions is not supported. New weight matrices in the objective function can be specified with the keyword arguments (see LinMPC
for the nomenclature). If Cwt ≠ Inf
, the augmented move suppression weight is $\mathbf{Ñ}_{H_c} = \mathrm{diag}(\mathbf{N}_{H_c}, C)$, else $\mathbf{Ñ}_{H_c} = \mathbf{N}_{H_c}$. The StateEstimator
mpc.estim
cannot be a Luenberger
observer or a SteadyKalmanFilter
(the default estimator). Construct the mpc
object with a time-varying KalmanFilter
instead. Note that the model is constant over the prediction horizon $H_p$.
Arguments
Keyword arguments with emphasis
are non-Unicode alternatives.
mpc::PredictiveController
: controller to set model and weights.model=mpc.estim.model
: new plant model (not supported byNonLinModel
).Mwt=nothing
: new main diagonal in $\mathbf{M}$ weight matrix (vector).Nwt=nothing
: new main diagonal in $\mathbf{N}$ weight matrix (vector).Lwt=nothing
: new main diagonal in $\mathbf{L}$ weight matrix (vector).M_Hp=nothing
: new $\mathbf{M}_{H_p}$ weight matrix.Ñ_Hc=nothing
orNtilde_Hc
: new $\mathbf{Ñ}_{H_c}$ weight matrix (see def. above).L_Hp=nothing
: new $\mathbf{L}_{H_p}$ weight matrix.- additional keyword arguments are passed to
setmodel!(mpc.estim)
.
Examples
julia> mpc = LinMPC(KalmanFilter(LinModel(ss(0.1, 0.5, 1, 0, 4.0)), σR=[√25]), Hp=1, Hc=1);
julia> mpc.estim.model.A[1], mpc.estim.R̂[1], mpc.weights.M_Hp[1], mpc.weights.Ñ_Hc[1]
(0.1, 25.0, 1.0, 0.1)
julia> setmodel!(mpc, LinModel(ss(0.42, 0.5, 1, 0, 4.0)); R̂=[9], M_Hp=[10], Nwt=[0.666]);
julia> mpc.estim.model.A[1], mpc.estim.R̂[1], mpc.weights.M_Hp[1], mpc.weights.Ñ_Hc[1]
(0.42, 9.0, 10.0, 0.666)
Get Additional Information
ModelPredictiveControl.getinfo
— Functiongetinfo(estim::MovingHorizonEstimator) -> info
Get additional info on estim
MovingHorizonEstimator
optimum for troubleshooting.
If estim.direct==true
, the function should be called after calling preparestate!
. Otherwise, call it after updatestate!
. It returns the dictionary info
with the following fields:
Fields with emphasis
are non-Unicode alternatives.
:Ŵ
or:What
: optimal estimated process noise over $N_k$, $\mathbf{Ŵ}$:ϵ
or:epsilon
: optimal slack variable, $ϵ$:X̂
or:Xhat
: optimal estimated states over $N_k+1$, $\mathbf{X̂}$:x̂
or:xhat
: optimal estimated state for the next time step, $\mathbf{x̂}_k(k+1)$:V̂
or:Vhat
: optimal estimated sensor noise over $N_k$, $\mathbf{V̂}$:P̄
or:Pbar
: estimation error covariance at arrival, $\mathbf{P̄}$:x̄
or:xbar
: optimal estimation error at arrival, $\mathbf{x̄}$:Ŷ
or:Yhat
: optimal estimated outputs over $N_k$, $\mathbf{Ŷ}$:Ŷm
or:Yhatm
: optimal estimated measured outputs over $N_k$, $\mathbf{Ŷ^m}$:x̂arr
or:xhatarr
: optimal estimated state at arrival, $\mathbf{x̂}_k(k-N_k+p)$:J
: objective value optimum, $J$:Ym
: measured outputs over $N_k$, $\mathbf{Y^m}$:U
: manipulated inputs over $N_k$, $\mathbf{U}$:D
: measured disturbances over $N_k$, $\mathbf{D}$:sol
: solution summary of the optimizer for printing
Examples
julia> model = LinModel(ss(1.0, 1.0, 1.0, 0, 5.0));
julia> estim = MovingHorizonEstimator(model, He=1, nint_ym=0, direct=false);
julia> updatestate!(estim, [0], [1]);
julia> round.(getinfo(estim)[:Ŷ], 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:
Fields with emphasis
are non-Unicode alternatives.
:ΔU
or:DeltaU
: optimal manipulated input increments over $H_c$, $\mathbf{ΔU}$:ϵ
or:epsilon
: optimal slack variable, $ϵ$:D̂
or:Dhat
: predicted measured disturbances over $H_p$, $\mathbf{D̂}$:ŷ
or:yhat
: current estimated output, $\mathbf{ŷ}(k)$:Ŷ
or:Yhat
: optimal predicted outputs over $H_p$, $\mathbf{Ŷ}$:Ŷs
or:Yhats
: predicted stochastic output over $H_p$ ofInternalModel
, $\mathbf{Ŷ_s}$:R̂y
or:Rhaty
: predicted output setpoint over $H_p$, $\mathbf{R̂_y}$:R̂u
or:Rhatu
: predicted manipulated input setpoint over $H_p$, $\mathbf{R̂_u}$:x̂end
or:xhatend
: optimal terminal states, $\mathbf{x̂}_i(k+H_p)$:J
: objective value optimum, $J$:U
: optimal manipulated inputs over $H_p$, $\mathbf{U}$:u
: current optimal manipulated input, $\mathbf{u}(k)$:d
: current measured disturbance, $\mathbf{d}(k)$
For LinMPC
and NonLinMPC
, the field :sol
also contains the optimizer solution summary that can be printed. Lastly, the economical cost :JE
and the custom nonlinear constraints :gc
values at the optimum are also available for NonLinMPC
.
Examples
julia> mpc = LinMPC(LinModel(tf(5, [2, 1]), 3), Nwt=[0], Hp=1, Hc=1);
julia> preparestate!(mpc, [0]); u = moveinput!(mpc, [10]);
julia> round.(getinfo(mpc)[:Ŷ], digits=3)
1-element Vector{Float64}:
10.0
Real-Time Simulate and Control
These utilities are for soft real-time applications. They are not suitable for hard real-time environnement like safety-critical processes.
Save current time t
ModelPredictiveControl.savetime!
— Functionsavetime!(model::SimModel) -> t
Set model.t
to time()
and return the value.
Used in conjunction with periodsleep
for simple soft real-time simulations. Call this function before any other in the simulation loop.
savetime!(estim::StateEstimator) -> t
Call savetime!(estim.model)
and return the time t
.
savetime!(mpc::PredictiveController) -> t
Call savetime!(mpc.estim.model)
and return the time t
.
Period Sleep
ModelPredictiveControl.periodsleep
— Functionperiodsleep(model::SimModel, busywait=false) -> nothing
Sleep for model.Ts
s minus the time elapsed since the last call to savetime!
.
It calls sleep
if busywait
is false
. Else, a simple while
loop implements busy-waiting. As a rule-of-thumb, busy-waiting should be used if model.Ts < 0.1
s, since the accuracy of sleep
is around 1 ms. Can be used to implement simple soft real-time simulations, see the example below.
The allocations in Julia are garbage-collected (GC) automatically. This can affect the timing. In such cases, you can temporarily stop the GC with GC.enable(false)
, and restart it at a convenient time e.g.: just before calling periodsleep
.
Examples
julia> model = LinModel(tf(2, [0.3, 1]), 0.25);
julia> function sim_realtime!(model)
t_0 = time()
for i=1:3
t = savetime!(model) # first function called
println(round(t - t_0, digits=3))
updatestate!(model, [1])
periodsleep(model, true) # last function called
end
end;
julia> sim_realtime!(model)
0.0
0.25
0.5
periodsleep(estim::StateEstimator) -> nothing
Call periodsleep(estim.model)
.
periodsleep(mpc::PredictiveController) -> nothing
Call periodsleep(mpc.estim.model)
.