Functions: Plant Models

The SimModel types represents discrete state-space models that can be used to construct StateEstimators and PredictiveControllers, or as plant simulators by calling evaloutput and updatestate! methods on SimModel instances (to test estimator/controller designs). For time simulations, the states $\mathbf{x}$ are stored inside SimModel instances. Use setstate! method to manually modify them.

SimModel

ModelPredictiveControl.SimModelType

Abstract supertype of LinModel and NonLinModel types.


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

Functor allowing callable SimModel object as an alias for evaloutput.

Examples

julia> model = NonLinModel((x,u,_,_)->-x + u, (x,_,_)->x .+ 20, 4, 1, 1, 1, solver=nothing);

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

LinModel

ModelPredictiveControl.LinModelType
LinModel(sys::StateSpace[, Ts]; i_u=1:size(sys,2), i_d=Int[])

Construct a linear model from state-space model sys with sampling time Ts in second.

The system sys can be continuous or discrete-time (Ts can be omitted for the latter). For continuous dynamics, its state-space equations are (discrete case in Extended Help):

\[\begin{aligned} \mathbf{ẋ}(t) &= \mathbf{A x}(t) + \mathbf{B z}(t) \\ \mathbf{y}(t) &= \mathbf{C x}(t) + \mathbf{D z}(t) \end{aligned}\]

with the state $\mathbf{x}$ and output $\mathbf{y}$ vectors. The $\mathbf{z}$ vector comprises the manipulated inputs $\mathbf{u}$ and measured disturbances $\mathbf{d}$, in any order. i_u provides the indices of $\mathbf{z}$ that are manipulated, and i_d, the measured disturbances. The constructor automatically discretizes continuous systems, resamples discrete ones if Ts ≠ sys.Ts, computes a new realization if not minimal, and separates the $\mathbf{z}$ terms in two parts (details in Extended Help). The rest of the documentation assumes discrete dynamics since all systems end up in this form.

See also ss

Examples

julia> model1 = LinModel(ss(-0.1, 1.0, 1.0, 0), 2.0) # continuous-time StateSpace
LinModel with a sample time Ts = 2.0 s and:
 1 manipulated inputs u
 1 states x
 1 outputs y
 0 measured disturbances d

julia> model2 = LinModel(ss(0.4, 0.2, 0.3, 0, 0.1)) # discrete-time StateSpace
LinModel with a sample time Ts = 0.1 s and:
 1 manipulated inputs u
 1 states x
 1 outputs y
 0 measured disturbances d

Extended Help

Extended Help

The state-space equations are similar if sys is discrete-time:

\[\begin{aligned} \mathbf{x}(k+1) &= \mathbf{A x}(k) + \mathbf{B z}(k) \\ \mathbf{y}(k) &= \mathbf{C x}(k) + \mathbf{D z}(k) \end{aligned}\]

Continuous dynamics are internally discretized using c2d and :zoh for manipulated inputs, and :tustin, for measured disturbances. Lastly, if sys is discrete and the provided argument Ts ≠ sys.Ts, the system is resampled by using the aforementioned discretization methods.

Note that the constructor transforms the system to its minimal realization using minreal for controllability/observability. As a consequence, the final state-space representation may be different from the one provided in sys. It is also converted into a more practical form ($\mathbf{D_u=0}$ because of the zero-order hold):

\[\begin{aligned} \mathbf{x}(k+1) &= \mathbf{A x}(k) + \mathbf{B_u u}(k) + \mathbf{B_d d}(k) \\ \mathbf{y}(k) &= \mathbf{C x}(k) + \mathbf{D_d d}(k) \end{aligned}\]

Use the syntax LinModel{NT}(A, Bu, C, Bd, Dd, Ts) to force a specific state-space representation.

source
LinModel(sys::TransferFunction[, Ts]; i_u=1:size(sys,2), i_d=Int[])

Convert to minimal realization state-space when sys is a transfer function.

sys is equal to $\frac{\mathbf{y}(s)}{\mathbf{z}(s)}$ for continuous-time, and $\frac{\mathbf{y}(z)}{\mathbf{z}(z)}$, for discrete-time.

See also tf

Examples

julia> model = LinModel([tf(3, [30, 1]) tf(-2, [5, 1])], 0.5, i_d=[2])
LinModel with a sample time Ts = 0.5 s and:
 1 manipulated inputs u
 2 states x
 1 outputs y
 1 measured disturbances d
source
LinModel(sys::DelayLtiSystem, Ts; i_u=1:size(sys,2), i_d=Int[])

Discretize with zero-order hold when sys is a continuous system with delays.

The delays must be multiples of the sample time Ts.

source
LinModel{NT}(A, Bu, C, Bd, Dd, Ts)

Construct the model from the discrete state-space matrices A, Bu, C, Bd, Dd directly.

See LinModel(::StateSpace) Extended Help for the meaning of the matrices. The arguments Bd and Dd can be the scalar 0 if there is no measured disturbance. This syntax do not modify the state-space representation provided in argument (minreal is not called). Care must be taken to ensure that the model is controllable and observable. The optional parameter NT explicitly set the number type of vectors (default to Float64).

source
LinModel(model::NonLinModel; x=model.x0+model.xop, u=model.uop, d=model.dop)

Call linearize(model; x, u, d) and return the resulting linear model.

source

NonLinModel

ModelPredictiveControl.NonLinModelType
NonLinModel{NT}(f::Function,  h::Function,  Ts, nu, nx, ny, nd=0; p=[], solver=RungeKutta(4))
NonLinModel{NT}(f!::Function, h!::Function, Ts, nu, nx, ny, nd=0; p=[], solver=RungeKutta(4))

Construct a nonlinear model from state-space functions f/f! and h/h!.

Both continuous and discrete-time models are supported. The default arguments assume continuous dynamics. Use solver=nothing for the discrete case (see Extended Help). The functions are defined as:

\[\begin{aligned} \mathbf{ẋ}(t) &= \mathbf{f}\Big( \mathbf{x}(t), \mathbf{u}(t), \mathbf{d}(t), \mathbf{p} \Big) \\ \mathbf{y}(t) &= \mathbf{h}\Big( \mathbf{x}(t), \mathbf{d}(t), \mathbf{p} \Big) \end{aligned}\]

where $\mathbf{x}$, $\mathbf{y}$, $\mathbf{u}$, $\mathbf{d}$ and $\mathbf{p}$ are respectively the state, output, manipulated input, measured disturbance and parameter vectors. If the dynamics is a function of the time, simply add a measured disturbance defined as $d(t) = t$. The functions can be implemented in two possible ways:

  1. Non-mutating functions (out-of-place): define them as f(x, u, d, p) -> ẋ and h(x, d, p) -> y. This syntax is simple and intuitive but it allocates more memory.
  2. Mutating functions (in-place): define them as f!(ẋ, x, u, d, p) -> nothing and h!(y, x, d, p) -> nothing. This syntax reduces the allocations and potentially the computational burden as well.

Ts is the sampling time in second. nu, nx, ny and nd are the respective number of manipulated inputs, states, outputs and measured disturbances. The keyword argument p is the parameters of the model passed to the two functions. It can be of any Julia object but use a mutable type if you want to change them later e.g.: a vector.

Tip

Replace the d or p argument with _ in your functions if not needed (see Examples below).

A 4th order RungeKutta solver discretizes the differential equations by default. The rest of the documentation assumes discrete dynamics since all models end up in this form. The optional parameter NT explicitly set the number type of vectors (default to Float64).

Warning

The two functions must be in pure Julia to use the model in NonLinMPC, ExtendedKalmanFilter, MovingHorizonEstimator and linearize.

See also LinModel.

Examples

julia> f!(ẋ, x, u, _ , p) = (ẋ .= p*x .+ u; nothing);

julia> h!(y, x, _ , _ ) = (y .= 0.1x; nothing);

julia> model1 = NonLinModel(f!, h!, 5.0, 1, 1, 1, p=-0.2)       # continuous dynamics
NonLinModel with a sample time Ts = 5.0 s, RungeKutta solver and:
 1 manipulated inputs u
 1 states x
 1 outputs y
 0 measured disturbances d

julia> f(x, u, _ , _ ) = 0.1x + u;

julia> h(x, _ , _ ) = 2x;

julia> model2 = NonLinModel(f, h, 2.0, 1, 1, 1, solver=nothing) # discrete dynamics
NonLinModel with a sample time Ts = 2.0 s, empty solver and:
 1 manipulated inputs u
 1 states x
 1 outputs y
 0 measured disturbances d

Extended Help

Extended Help

State-space functions are similar for discrete dynamics:

\[\begin{aligned} \mathbf{x}(k+1) &= \mathbf{f}\Big( \mathbf{x}(k), \mathbf{u}(k), \mathbf{d}(k), \mathbf{p} \Big) \\ \mathbf{y}(k) &= \mathbf{h}\Big( \mathbf{x}(k), \mathbf{d}(k), \mathbf{p} \Big) \end{aligned}\]

with two possible implementations as well:

  1. Non-mutating functions: define them as f(x, u, d, p) -> xnext and h(x, d, p) -> y.
  2. Mutating functions: define them as f!(xnext, x, u, d, p) -> nothing and h!(y, x, d, p) -> nothing.
source

Set Variable Names

ModelPredictiveControl.setname!Function
setname!(model::SimModel; u=nothing, y=nothing, d=nothing, x=nothing) -> model

Set the names of model inputs u, outputs y, disturbances d, and states x.

The keyword arguments u, y, d, and x must be vectors of strings. The strings are used in the plotting functions.

Examples

julia> model = setname!(LinModel(tf(3, [10, 1]), 2.0), u=["\$A\$ (%)"], y=["\$T\$ (∘C)"])
LinModel with a sample time Ts = 2.0 s and:
 1 manipulated inputs u
 1 states x
 1 outputs y
 0 measured disturbances d
source

Set Operating Points

ModelPredictiveControl.setop!Function
setop!(model; uop=nothing, yop=nothing, dop=nothing, xop=nothing, fop=nothing) -> model

Set the operating points of model (both LinModel and NonLinModel).

Introducing deviations vectors around manipulated input uop, model output yop, measured disturbance dop, and model state xop operating points (a.k.a. nominal values):

\[\begin{align*} \mathbf{u_0}(k) &= \mathbf{u}(k) - \mathbf{u_{op}} \\ \mathbf{d_0}(k) &= \mathbf{d}(k) - \mathbf{d_{op}} \\ \mathbf{y_0}(k) &= \mathbf{y}(k) - \mathbf{y_{op}} \\ \mathbf{x_0}(k) &= \mathbf{x}(k) - \mathbf{x_{op}} \\ \end{align*}\]

The state-space description of LinModel around the operating points is:

\[\begin{aligned} \mathbf{x_0}(k+1) &= \mathbf{A x_0}(k) + \mathbf{B_u u_0}(k) + \mathbf{B_d d_0}(k) + \mathbf{f_{op}} - \mathbf{x_{op}} \\ \mathbf{y_0}(k) &= \mathbf{C x_0}(k) + \mathbf{D_d d_0}(k) \end{aligned}\]

and, for NonLinModel:

\[\begin{aligned} \mathbf{x_0}(k+1) &= \mathbf{f}\Big(\mathbf{x_0}(k), \mathbf{u_0}(k), \mathbf{d_0}(k), \mathbf{p}\Big) + \mathbf{f_{op}} - \mathbf{x_{op}} \\ \mathbf{y_0}(k) &= \mathbf{h}\Big(\mathbf{x_0}(k), \mathbf{d_0}(k), \mathbf{p}\Big) \end{aligned}\]

The state xop and the additional fop operating points are frequently zero e.g.: when model comes from system identification. The two vectors are internally used by linearize for non-equilibrium points.

Examples

julia> model = setop!(LinModel(tf(3, [10, 1]), 2.0), uop=[50], yop=[20])
LinModel with a sample time Ts = 2.0 s and:
 1 manipulated inputs u
 1 states x
 1 outputs y
 0 measured disturbances d

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

Linearize

ModelPredictiveControl.linearizeFunction
linearize(model::SimModel; x=model.x0+model.xop, u=model.uop, d=model.dop) -> linmodel

Linearize model at the operating points x, u, d and return the LinModel.

The arguments x, u and d are the linearization points for the state $\mathbf{x}$, manipulated input $\mathbf{u}$ and measured disturbance $\mathbf{d}$, respectively (not necessarily an equilibrium, details in Extended Help). The Jacobians of $\mathbf{f}$ and $\mathbf{h}$ functions are automatically computed with ForwardDiff.jl.

Warning

See Extended Help if you get an error like: MethodError: no method matching (::var"##")(::Vector{ForwardDiff.Dual}).

Examples

julia> model = NonLinModel((x,u,_,_)->x.^3 + u, (x,_,_)->x, 0.1, 1, 1, 1, solver=nothing);

julia> linmodel = linearize(model, x=[10.0], u=[0.0]); 

julia> linmodel.A
1×1 Matrix{Float64}:
 300.0

Extended Help

Extended Help

With the nonlinear state-space model:

\[\begin{aligned} \mathbf{x}(k+1) &= \mathbf{f}\Big(\mathbf{x}(k), \mathbf{u}(k), \mathbf{d}(k), \mathbf{p}\Big) \\ \mathbf{y}(k) &= \mathbf{h}\Big(\mathbf{x}(k), \mathbf{d}(k), \mathbf{p}\Big) \end{aligned}\]

its linearization at the operating point $\mathbf{x_{op}, u_{op}, d_{op}}$ is:

\[\begin{aligned} \mathbf{x_0}(k+1) &≈ \mathbf{A x_0}(k) + \mathbf{B_u u_0}(k) + \mathbf{B_d d_0}(k) + \mathbf{f(x_{op}, u_{op}, d_{op}, p)} - \mathbf{x_{op}} \\ \mathbf{y_0}(k) &≈ \mathbf{C x_0}(k) + \mathbf{D_d d_0}(k) \end{aligned}\]

based on the deviation vectors $\mathbf{x_0, u_0, d_0, y_0}$ introduced in setop! documentation, and the Jacobian matrices:

\[\begin{aligned} \mathbf{A} &= \left. \frac{∂\mathbf{f(x, u, d, p)}}{∂\mathbf{x}} \right|_{\mathbf{x=x_{op},\, u=u_{op},\, d=d_{op}}} \\ \mathbf{B_u} &= \left. \frac{∂\mathbf{f(x, u, d, p)}}{∂\mathbf{u}} \right|_{\mathbf{x=x_{op},\, u=u_{op},\, d=d_{op}}} \\ \mathbf{B_d} &= \left. \frac{∂\mathbf{f(x, u, d, p)}}{∂\mathbf{d}} \right|_{\mathbf{x=x_{op},\, u=u_{op},\, d=d_{op}}} \\ \mathbf{C} &= \left. \frac{∂\mathbf{h(x, d, p)}}{∂\mathbf{x}} \right|_{\mathbf{x=x_{op},\, d=d_{op}}} \\ \mathbf{D_d} &= \left. \frac{∂\mathbf{h(x, d, p)}}{∂\mathbf{d}} \right|_{\mathbf{x=x_{op},\, d=d_{op}}} \end{aligned}\]

Following setop! notation, we find:

\[\begin{aligned} \mathbf{f_{op}} &= \mathbf{f(x_{op}, u_{op}, d_{op}, p)} \\ \mathbf{y_{op}} &= \mathbf{h(x_{op}, d_{op}, p)} \end{aligned}\]

Notice that $\mathbf{f_{op} - x_{op} = 0}$ if the point is an equilibrium. The equations are similar if the nonlinear model has nonzero operating points.

Automatic differentiation (AD) allows exact Jacobians. The NonLinModel f and h functions must be compatible with this feature though. See Automatic differentiation for common mistakes when writing these functions.

source
ModelPredictiveControl.linearize!Function
linearize!(linmodel::LinModel, model::SimModel; <keyword arguments>) -> linmodel

Linearize model and store the result in linmodel (in-place).

The keyword arguments are identical to linearize.

Examples

julia> model = NonLinModel((x,u,_,_)->x.^3 + u, (x,_,_)->x, 0.1, 1, 1, 1, solver=nothing);

julia> linmodel = linearize(model, x=[10.0], u=[0.0]); linmodel.A
1×1 Matrix{Float64}:
 300.0

julia> linearize!(linmodel, model, x=[20.0], u=[0.0]); linmodel.A
1×1 Matrix{Float64}:
 1200.0
source

Differential Equation Solvers

DiffSolver

RungeKutta

ModelPredictiveControl.RungeKuttaType
RungeKutta(order=4; supersample=1)

Create a Runge-Kutta solver with optional super-sampling.

Only the 4th order Runge-Kutta is supported for now. The keyword argument supersample provides the number of internal steps (default to 1 step).

source