doc: modify NonLinMPC pendulum example in manual with new p argument

Author: franckgaga

docs/src/manual/nonlinmpc.md

@@ -40,34 +40,34 @@ The plant model is nonlinear:
 
 in which ``g`` is the gravitational acceleration in m/s², ``L``, the pendulum length in m,
 ``K``, the friction coefficient at the pivot point in kg/s, and ``m``, the mass attached at
-the end of the pendulum in kg. The [`NonLinModel`](@ref) constructor assumes by default
-that the state function `f` is continuous in time, that is, an ordinary differential
-equation system (like here):
+the end of the pendulum in kg, all bundled in the parameter vector ``\mathbf{p} =
+[\begin{smallmatrix} g & L & K & m \end{smallmatrix}]'``. The [`NonLinModel`](@ref)
+constructor assumes by default that the state function `f` is continuous in time, that is,
+an ordinary differential equation system (like here):
 
 ```@example 1
 using ModelPredictiveControl
-function pendulum(par, x, u)
-    g, L, K, m = par        # [m/s²], [m], [kg/s], [kg]
+function f(x, u, _ , p)
+    g, L, K, m = p          # [m/s²], [m], [kg/s], [kg]
     θ, ω = x[1], x[2]       # [rad], [rad/s]
     τ  = u[1]               # [Nm]
     dθ = ω
     dω = -g/L*sin(θ) - K/m*ω + τ/m/L^2
     return [dθ, dω]
 end
-# declared constants, to avoid type-instability in the f function, for speed:
-const par = (9.8, 0.4, 1.2, 0.3)
-f(x, u, _ ) = pendulum(par, x, u)
-h(x, _ )    = [180/π*x[1]]  # [°]
+h(x, _ , _ ) = [180/π*x[1]] # [°]
+p_model = [9.8, 0.4, 1.2, 0.3]
 nu, nx, ny, Ts = 1, 2, 1, 0.1
 vu, vx, vy = ["\$τ\$ (Nm)"], ["\$θ\$ (rad)", "\$ω\$ (rad/s)"], ["\$θ\$ (°)"]
-model = setname!(NonLinModel(f, h, Ts, nu, nx, ny); u=vu, x=vx, y=vy)
+model = setname!(NonLinModel(f, h, Ts, nu, nx, ny; p=p_model); u=vu, x=vx, y=vy)
 ```
 
 The output function ``\mathbf{h}`` converts the ``θ`` angle to degrees. Note that special
 characters like ``θ`` can be typed in the Julia REPL or VS Code by typing `\theta` and
-pressing the `<TAB>` key. The tuple `par` is constant here to improve the [performance](https://docs.julialang.org/en/v1/manual/performance-tips/#Avoid-untyped-global-variables).
-A 4th order [`RungeKutta`](@ref) method solves the differential equations by default. It is
-good practice to first simulate `model` using [`sim!`](@ref) as a quick sanity check:
+pressing the `<TAB>` key. Note that the model parameter `p` can be any 
+ A 4th order [`RungeKutta`](@ref) method solves the differential
+equations by default. It is good practice to first simulate `model` using [`sim!`](@ref) as
+a quick sanity check:
 
 ```@example 1
 using Plots
@@ -101,9 +101,9 @@ motor torque ``τ``, with an associated standard deviation `σQint_u` of 0.1 N m
 estimator tuning is tested on a plant with a 25 % larger friction coefficient ``K``:
 
 ```@example 1
-const par_plant = (par[1], par[2], 1.25*par[3], par[4])
-f_plant(x, u, _ ) = pendulum(par_plant, x, u)
-plant = setname!(NonLinModel(f_plant, h, Ts, nu, nx, ny); u=vu, x=vx, y=vy)
+p_plant = copy(p_model)
+p_plant[3] = 1.25*p_model[3]
+plant = setname!(NonLinModel(f, h, Ts, nu, nx, ny; p=p_plant); u=vu, x=vx, y=vy)
 res = sim!(estim, N, [0.5], plant=plant, y_noise=[0.5])
 plot(res, plotu=false, plotxwithx̂=true)
 savefig("plot2_NonLinMPC.svg"); nothing # hide
@@ -182,10 +182,10 @@ angle ``θ`` is measured here). As the arguments of [`NonLinMPC`](@ref) economic
 Kalman Filter similar to the previous one (``\mathbf{y^m} = θ`` and ``\mathbf{y^u} = ω``):
 
 ```@example 1
-h2(x, _ ) = [180/π*x[1], x[2]]
+h2(x, _ , _ ) = [180/π*x[1], x[2]]
 nu, nx, ny = 1, 2, 2
-model2 = setname!(NonLinModel(f      , h2, Ts, nu, nx, ny), u=vu, x=vx, y=[vy; vx[2]])
-plant2 = setname!(NonLinModel(f_plant, h2, Ts, nu, nx, ny), u=vu, x=vx, y=[vy; vx[2]])
+model2 = setname!(NonLinModel(f, h2, Ts, nu, nx, ny; p=p_model), u=vu, x=vx, y=[vy; vx[2]])
+plant2 = setname!(NonLinModel(f, h2, Ts, nu, nx, ny; p=p_plant), u=vu, x=vx, y=[vy; vx[2]])
 estim2 = UnscentedKalmanFilter(model2; σQ, σR, nint_u, σQint_u, i_ym=[1])
 ```
 

src/controller/nonlinmpc.jl

@@ -159,7 +159,7 @@ This method uses the default state estimator :
 
 # Examples
 ```jldoctest
-julia> model = NonLinModel((x,u,_)->0.5x+u, (x,_)->2x, 10.0, 1, 1, 1, solver=nothing);
+julia> model = NonLinModel((x,u,_,_)->0.5x+u, (x,_,_)->2x, 10.0, 1, 1, 1, solver=nothing);
 
 julia> mpc = NonLinMPC(model, Hp=20, Hc=1, Cwt=1e6)
 NonLinMPC controller with a sample time Ts = 10.0 s, Ipopt optimizer, UnscentedKalmanFilter estimator and:
@@ -236,7 +236,7 @@ Use custom state estimator `estim` to construct `NonLinMPC`.
 
 # Examples
 ```jldoctest
-julia> model = NonLinModel((x,u,_)->0.5x+u, (x,_)->2x, 10.0, 1, 1, 1, solver=nothing);
+julia> model = NonLinModel((x,u,_,_)->0.5x+u, (x,_,_)->2x, 10.0, 1, 1, 1, solver=nothing);
 
 julia> estim = UnscentedKalmanFilter(model, σQint_ym=[0.05]);
 

src/estimator/kalman.jl

@@ -623,7 +623,7 @@ initial state estimate ``\mathbf{x̂}_{-1}(0)`` can be manually specified with [
 
 # Examples
 ```jldoctest
-julia> model = NonLinModel((x,u,_)->0.1x+u, (x,_)->2x, 10.0, 1, 1, 1, solver=nothing);
+julia> model = NonLinModel((x,u,_,_)->0.1x+u, (x,_,_)->2x, 10.0, 1, 1, 1, solver=nothing);
 
 julia> estim = UnscentedKalmanFilter(model, σR=[1], nint_ym=[2], σPint_ym_0=[1, 1])
 UnscentedKalmanFilter estimator with a sample time Ts = 10.0 s, NonLinModel and:
@@ -951,7 +951,7 @@ automatic differentiation.
 
 # Examples
 ```jldoctest
-julia> model = NonLinModel((x,u,_)->0.2x+u, (x,_)->-3x, 5.0, 1, 1, 1, solver=nothing);
+julia> model = NonLinModel((x,u,_,_)->0.2x+u, (x,_,_)->-3x, 5.0, 1, 1, 1, solver=nothing);
 
 julia> estim = ExtendedKalmanFilter(model, σQ=[2], σQint_ym=[2], σP_0=[0.1], σPint_ym_0=[0.1])
 ExtendedKalmanFilter estimator with a sample time Ts = 5.0 s, NonLinModel and:

src/estimator/mhe/construct.jl

@@ -256,7 +256,7 @@ See [`UnscentedKalmanFilter`](@ref) for details on the augmented process model a
 
 # Examples
 ```jldoctest
-julia> model = NonLinModel((x,u,_)->0.1x+u, (x,_)->2x, 10.0, 1, 1, 1, solver=nothing);
+julia> model = NonLinModel((x,u,_,_)->0.1x+u, (x,_,_)->2x, 10.0, 1, 1, 1, solver=nothing);
 
 julia> estim = MovingHorizonEstimator(model, He=5, σR=[1], σP_0=[0.01])
 MovingHorizonEstimator estimator with a sample time Ts = 10.0 s, Ipopt optimizer, NonLinModel and:

src/model/linearization.jl

@@ -21,7 +21,7 @@ necessarily an equilibrium, details in Extended Help). The Jacobians of ``\mathb
 
 # Examples
 ```jldoctest
-julia> model = NonLinModel((x,u,_)->x.^3 + u, (x,_)->x, 0.1, 1, 1, 1, solver=nothing);
+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]); 
 
@@ -96,7 +96,7 @@ The keyword arguments are identical to [`linearize`](@ref).
 
 # Examples
 ```jldoctest
-julia> model = NonLinModel((x,u,_)->x.^3 + u, (x,_)->x, 0.1, 1, 1, 1, solver=nothing);
+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}:

src/model/nonlinmodel.jl

@@ -69,7 +69,8 @@ functions are defined as:
 ```
 where ``\mathbf{x}``, ``\mathbf{y}, ``\mathbf{u}``, ``\mathbf{d}`` and ``\mathbf{p}`` are
 respectively the state, output, manipulated input, measured disturbance and parameter vectors,
-and ``t`` the time in second. The functions can be implemented in two possible ways:
+and ``t`` the time in second. If the dynamics is a function of 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.

src/plot_sim.jl

@@ -118,7 +118,7 @@ internal states.
 
 # Examples
 ```jldoctest
-julia> plant = NonLinModel((x,u,d)->0.1x+u+d, (x,_)->2x, 10.0, 1, 1, 1, 1, solver=nothing);
+julia> plant = NonLinModel((x,u,d,_)->0.1x+u+d, (x,_,_)->2x, 5, 1, 1, 1, 1, solver=nothing);
 
 julia> res = sim!(plant, 15, [0], [0], x_0=[1])
 Simulation results of NonLinModel with 15 time steps.

src/precompile.jl

@@ -74,8 +74,8 @@ exmpc.estim()
 u = exmpc([55, 30])
 sim!(exmpc, 3, [55, 30])
 
-f(x,u,_) = model.A*x + model.Bu*u
-h(x,_) = model.C*x
+f(x,u,_,_) = model.A*x + model.Bu*u
+h(x,_,_) = model.C*x
 
 nlmodel = setop!(NonLinModel(f, h, Ts, 2, 2, 2, solver=nothing), uop=[10, 10], yop=[50, 30])
 y = nlmodel()

src/sim_model.jl

@@ -11,7 +11,7 @@ Functor allowing callable `SimModel` object as an alias for [`evaloutput`](@ref)
 
 # Examples
 ```jldoctest
-julia> model = NonLinModel((x,u,_)->-x + u, (x,_)->x .+ 20, 10.0, 1, 1, 1, solver=nothing);
+julia> model = NonLinModel((x,u,_,_)->-x + u, (x,_,_)->x .+ 20, 4, 1, 1, 1, solver=nothing);
 
 julia> y = model()
 1-element Vector{Float64}:

test/test_predictive_control.jl

@@ -474,8 +474,8 @@ end
     linmodel1 = LinModel(sys,Ts,i_d=[3])
     nmpc0 = NonLinMPC(linmodel1, Hp=15)
     @test isa(nmpc0.estim, SteadyKalmanFilter)
-    f = (x,u,d) -> linmodel1.A*x + linmodel1.Bu*u + linmodel1.Bd*d
-    h = (x,d)   -> linmodel1.C*x + linmodel1.Dd*d
+    f = (x,u,d,_) -> linmodel1.A*x + linmodel1.Bu*u + linmodel1.Bd*d
+    h = (x,d,_)   -> linmodel1.C*x + linmodel1.Dd*d
     nonlinmodel = NonLinModel(f, h, Ts, 2, 4, 2, 1, solver=nothing)
     nmpc1 = NonLinMPC(nonlinmodel, Hp=15)
     @test isa(nmpc1.estim, UnscentedKalmanFilter)
@@ -545,8 +545,8 @@ end
     # ensure that the current estimated output is updated for correct JE values:
     @test nmpc.ŷ ≈ evaloutput(nmpc.estim, Float64[])
     linmodel2 = LinModel([tf(5, [2000, 1]) tf(7, [8000,1])], 3000.0, i_d=[2])
-    f = (x,u,d) -> linmodel2.A*x + linmodel2.Bu*u + linmodel2.Bd*d
-    h = (x,d)   -> linmodel2.C*x + linmodel2.Dd*d
+    f = (x,u,d,_) -> linmodel2.A*x + linmodel2.Bu*u + linmodel2.Bd*d
+    h = (x,d,_)   -> linmodel2.C*x + linmodel2.Dd*d
     nonlinmodel = NonLinModel(f, h, 3000.0, 1, 2, 1, 1, solver=nothing)
     nmpc2 = NonLinMPC(nonlinmodel, Nwt=[0], Hp=1000, Hc=1)
     d = [0.1]
@@ -628,8 +628,8 @@ end
 
 @testset "NonLinMPC other methods" begin
     linmodel = setop!(LinModel(sys,Ts,i_u=[1,2]), uop=[10,50], yop=[50,30])
-    f = (x,u,_) -> linmodel.A*x + linmodel.Bu*u
-    h = (x,_)   -> linmodel.C*x
+    f = (x,u,_,_) -> linmodel.A*x + linmodel.Bu*u
+    h = (x,_,_)   -> linmodel.C*x
     nonlinmodel = setop!(
         NonLinModel(f, h, Ts, 2, 2, 2, solver=nothing), uop=[10,50], yop=[50,30]
     )
@@ -652,8 +652,8 @@ end
     setconstraint!(nmpc_lin, c_ymin=[1.0,1.1], c_ymax=[1.2,1.3])
     @test all((-nmpc_lin.con.A_Ymin[:, end], -nmpc_lin.con.A_Ymax[:, end]) .≈ ([1.0,1.1], [1.2,1.3]))
 
-    f = (x,u,d) -> linmodel1.A*x + linmodel1.Bu*u + linmodel1.Bd*d
-    h = (x,d)   -> linmodel1.C*x + linmodel1.Dd*d
+    f = (x,u,d,_) -> linmodel1.A*x + linmodel1.Bu*u + linmodel1.Bd*d
+    h = (x,d,_)   -> linmodel1.C*x + linmodel1.Dd*d
     nonlinmodel = NonLinModel(f, h, Ts, 2, 4, 2, 1, solver=nothing)
     nmpc = NonLinMPC(nonlinmodel, Hp=1, Hc=1)
 
@@ -715,8 +715,8 @@ end
     info = getinfo(nmpc_lin)
     @test info[:x̂end][1] ≈ 0 atol=1e-1
 
-    f = (x,u,_) -> linmodel.A*x + linmodel.Bu*u
-    h = (x,_)   -> linmodel.C*x
+    f = (x,u,_,_) -> linmodel.A*x + linmodel.Bu*u
+    h = (x,_,_)   -> linmodel.C*x
     nonlinmodel = NonLinModel(f, h, linmodel.Ts, 1, 1, 1, solver=nothing)
     nmpc = NonLinMPC(nonlinmodel, Hp=50, Hc=5)
 
@@ -791,8 +791,8 @@ end
     @test mpc.M_Hp ≈ diagm(1:1000)
     @test mpc.Ñ_Hc ≈ diagm([0.1;1e6])
     @test mpc.L_Hp ≈ diagm(1.1:1000.1)
-    f = (x,u,d) -> estim.model.A*x + estim.model.Bu*u + estim.model.Bd*d
-    h = (x,d)   -> estim.model.C*x + estim.model.Du*d
+    f = (x,u,d,_) -> estim.model.A*x + estim.model.Bu*u + estim.model.Bd*d
+    h = (x,d,_)   -> estim.model.C*x + estim.model.Du*d
     nonlinmodel = NonLinModel(f, h, 10.0, 1, 1, 1)
     nmpc = NonLinMPC(nonlinmodel, Nwt=[0], Cwt=1e4, Hp=1000, Hc=10)
     setmodel!(nmpc, Mwt=[100], Nwt=[200], Lwt=[300])
@@ -808,8 +808,8 @@ end
 
 @testset "LinMPC v.s. NonLinMPC" begin
     linmodel = setop!(LinModel(sys,Ts,i_d=[3]), uop=[10,50], yop=[50,30], dop=[20])
-    f = (x,u,d) -> linmodel.A*x + linmodel.Bu*u + linmodel.Bd*d
-    h = (x,d)   -> linmodel.C*x + linmodel.Dd*d
+    f = (x,u,d,_) -> linmodel.A*x + linmodel.Bu*u + linmodel.Bd*d
+    h = (x,d,_)   -> linmodel.C*x + linmodel.Dd*d
     nonlinmodel = NonLinModel(f, h, Ts, 2, 4, 2, 1, solver=nothing)
     nonlinmodel = setop!(nonlinmodel, uop=[10,50], yop=[50,30], dop=[20])
     optim = JuMP.Model(optimizer_with_attributes(Ipopt.Optimizer, "sb"=>"yes"))