added : DAQP solver in LinMPC manual

Author: franckgaga

docs/Project.toml

@@ -1,9 +1,10 @@
 [deps]
-Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
 ControlSystemsBase = "aaaaaaaa-a6ca-5380-bf3e-84a91bcd477e"
+DAQP = "c47d62df-3981-49c8-9651-128b1cd08617"
+Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
 
 [compat]
 ControlSystemsBase = "1"
-Documenter = "0.27"
+Documenter = "0.27"

docs/src/manual/linmpc.md

@@ -68,11 +68,12 @@ We design our [`LinMPC`](@ref) controllers by including the level constraint wit
 mpc = setconstraint!(LinMPC(model, Hp=15, Hc=2), ŷmin=[45, -Inf])
 ```
 
-By default, [`LinMPC`](@ref) controllers use a [`SteadyKalmanFilter`](@ref) to estimate the
-plant states. Before closing the loop, we call [`initstate!`](@ref) with the actual plant
-inputs and measurements to ensure a bumpless transfer. Since `model` simulates our plant
-here, its output will initialize the states. [`LinModel`](@ref) objects are callable for
-this purpose (an alias for [`evaloutput`](@ref)):
+By default, [`LinMPC`](@ref) controllers use [`OSQP`](https://osqp.org/) to solve the
+problem and a [`SteadyKalmanFilter`](@ref) to estimate the plant states. Before closing the
+loop, we call [`initstate!`](@ref) with the actual plant inputs and measurements to ensure a
+bumpless transfer. Since `model` simulates our plant here, its output will initialize the
+states. [`LinModel`](@ref) objects are callable for this purpose (an alias for
+[`evaloutput`](@ref)):
 
 ```@example 1
 u = model.uop
@@ -120,12 +121,34 @@ Lastly, we plot the closed-loop test with the `Plots` package:
 
 ```@example 1
 using Plots
-p1 = plot(t_data, y_data[1,:], label="level"); ylabel!("level")
-plot!(t_data, ry_data[1,:], label="setpoint", linestyle=:dash, linetype=:steppost)
-plot!(t_data, fill(45,size(t_data)), label="min", linestyle=:dot, linetype=:steppost)
-p2 = plot(t_data, y_data[2,:], label="temp."); ylabel!("temp.")
-plot!(t_data, ry_data[2,:],label="setpoint", linestyle=:dash, linetype=:steppost)
-p3 = plot(t_data,u_data[1,:],label="cold", linetype=:steppost); ylabel!("flow rate")
-plot!(t_data,u_data[2,:],label="hot", linetype=:steppost); xlabel!("time (s)")
-p = plot(p1, p2, p3, layout=(3,1), fmt=:svg)
+function plot_data(t_data, u_data, y_data, ry_data)
+    p1 = plot(t_data, y_data[1,:], label="level"); ylabel!("level")
+    plot!(t_data, ry_data[1,:], label="setpoint", linestyle=:dash, linetype=:steppost)
+    plot!(t_data, fill(45,size(t_data)), label="min", linestyle=:dot, linewidth=1.5)
+    p2 = plot(t_data, y_data[2,:], label="temp."); ylabel!("temp.")
+    plot!(t_data, ry_data[2,:],label="setpoint", linestyle=:dash, linetype=:steppost)
+    p3 = plot(t_data,u_data[1,:],label="cold", linetype=:steppost); ylabel!("flow rate")
+    plot!(t_data,u_data[2,:],label="hot", linetype=:steppost); xlabel!("time (s)")
+    return plot(p1, p2, p3, layout=(3,1), fmt=:svg)
+end
+plot_data(t_data, u_data, y_data, ry_data)
+```
+
+For some situations, when [`LinMPC`](@ref) matrices are small/medium and dense, [`DAQP`](https://darnstrom.github.io/daqp/)
+optimizer may be more efficient. To install it, run:
+
+```julia
+using Pkg; Pkg.add("DAQP")
+```
+
+Constructing a [`LinMPC`](@ref) with `DAQP` leads to identical results here:
+
+```@example 1
+using JuMP, DAQP
+daqp = Model(DAQP.Optimizer)
+mpc2 = setconstraint!(LinMPC(model, Hp=15, Hc=2, optim=daqp), ŷmin=[45, -Inf])
+setstate!(model, zeros(model.nx))
+initstate!(mpc2, model.uop, model())
+u_data2, y_data2, ry_data2 = test_mpc(mpc2, model)
+plot_data(t_data, u_data2, y_data2, ry_data2)
 ```

example/juMPC.jl

@@ -100,7 +100,7 @@ setconstraint!(nmpc2, Δumin=[-Inf,-Inf],Δumax=[+Inf,+Inf])
 nx = linModel4.nx
 kf = KalmanFilter(linModel4, σP0=10*ones(nx), σQ=0.01*ones(nx), σR=[0.1, 0.1], σQ_int=0.05*ones(2), σP0_int=10*ones(2))
 
-mpc = LinMPC(kf, Hp=15, Hc=1, Mwt=[1, 1] , Nwt=[0.1, 0.1], Cwt=1e5)#, optim=Model(DAQP.Optimizer))#, optim=Model(HiGHS.Optimizer))#, optim=Model(DAQP.Optimizer))
+mpc = LinMPC(kf, Hp=15, Hc=1, Mwt=[1, 1] , Nwt=[0.1, 0.1], Cwt=1e5)#, optim=Model(HiGHS.Optimizer))#, optim=Model(DAQP.Optimizer))
 
 setconstraint!(mpc, c_umin=[0,0], c_umax=[0,0])
 setconstraint!(mpc, c_ŷmin=[1,1], c_ŷmax=[1,1])
@@ -164,9 +164,11 @@ using PlotThemes, Plots
 theme(:dark)
 default(fontfamily="Computer Modern"); scalefontsizes(1.1)
 
-test_mpc(linModel4 , mpc)
-@time u_data, y_data, r_data, d_data = test_mpc(linModel4, mpc)
+using BenchmarkTools
 
+test_mpc(linModel4 , mpc)
+@btime u_data, y_data, r_data, d_data = test_mpc(linModel4, mpc)
+#=
 resM = sim!(nonLinModel2, mpc.Hp+10, [1,-1])
 psM  = plot(resM)
 display(psM)
@@ -205,3 +207,4 @@ display(pd)
 display(pu)
 display(py)
 
+=#

src/controller/linmpc.jl

@@ -134,7 +134,7 @@ arguments.
 julia> model = LinModel([tf(3, [30, 1]); tf(-2, [5, 1])], 4);
 
 julia> mpc = LinMPC(model, Mwt=[0, 1], Nwt=[0.5], Hp=30, Hc=1)
-LinMPC controller with a sample time Ts = 4.0 s, SteadyKalmanFilter estimator and:
+LinMPC controller with a sample time Ts = 4.0 s, OSQP optimizer, SteadyKalmanFilter estimator and:
  30 prediction steps Hp
   1 control steps Hc
   1 manipulated inputs u
@@ -164,7 +164,7 @@ Use custom state estimator `estim` to construct `LinMPC`.
 julia> estim = KalmanFilter(LinModel([tf(3, [30, 1]); tf(-2, [5, 1])], 4), i_ym=[2]);
 
 julia> mpc = LinMPC(estim, Mwt=[0, 1], Nwt=[0.5], Hp=30, Hc=1)
-LinMPC controller with a sample time Ts = 4.0 s, KalmanFilter estimator and:
+LinMPC controller with a sample time Ts = 4.0 s, OSQP optimizer, KalmanFilter estimator and:
  30 prediction steps Hp
   1 control steps Hc
   1 manipulated inputs u

src/controller/nonlinmpc.jl

@@ -134,7 +134,7 @@ This method uses the default state estimator :
 julia> model = NonLinModel((x,u,_)->0.5x+u, (x,_)->2x, 10.0, 1, 1, 1);
 
 julia> mpc = NonLinMPC(model, Hp=20, Hc=1, Cwt=1e6)
-NonLinMPC controller with a sample time Ts = 10.0 s, UnscentedKalmanFilter estimator and:
+NonLinMPC controller with a sample time Ts = 10.0 s, Ipopt optimizer, UnscentedKalmanFilter estimator and:
  20 prediction steps Hp
   1 control steps Hc
   1 manipulated inputs u
@@ -171,7 +171,7 @@ julia> model = NonLinModel((x,u,_)->0.5x+u, (x,_)->2x, 10.0, 1, 1, 1);
 julia> estim = UnscentedKalmanFilter(model, σQ_int=[0.05]);
 
 julia> mpc = NonLinMPC(estim, Hp=20, Hc=1, Cwt=1e6)
-NonLinMPC controller with a sample time Ts = 10.0 s, UnscentedKalmanFilter estimator and:
+NonLinMPC controller with a sample time Ts = 10.0 s, Ipopt optimizer, UnscentedKalmanFilter estimator and:
  20 prediction steps Hp
   1 control steps Hc
   1 manipulated inputs u

src/predictive_control.jl

@@ -93,7 +93,7 @@ julia> mpc = LinMPC(setop!(LinModel(tf(3, [30, 1]), 4), uop=[50], yop=[25]));
 julia> mpc = setconstraint!(mpc, umin=[0], umax=[100], c_umin=[0.0], c_umax=[0.0]);
 
 julia> mpc = setconstraint!(mpc, Δumin=[-10], Δumax=[+10], c_Δumin=[1.0], c_Δumax=[1.0])
-LinMPC controller with a sample time Ts = 4.0 s, SteadyKalmanFilter estimator and:
+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
@@ -1009,7 +1009,8 @@ function Base.show(io::IO, mpc::PredictiveController)
     nx̂, nym, nyu = mpc.estim.nx̂, mpc.estim.nym, mpc.estim.nyu
     n = maximum(ndigits.((Hp, Hc, nu, nx̂, nym, nyu, nd))) + 1
     println(io, "$(typeof(mpc).name.name) controller with a sample time Ts = "*
-                "$(mpc.estim.model.Ts) s, $(typeof(mpc.estim).name.name) estimator and:")
+                "$(mpc.estim.model.Ts) s, $(solver_name(mpc.optim)) optimizer, "*
+                "$(typeof(mpc.estim).name.name) estimator and:")
     println(io, "$(lpad(Hp, n)) prediction steps Hp")
     println(io, "$(lpad(Hc, n)) control steps Hc")
     println(io, "$(lpad(nu, n)) manipulated inputs u")