diff --git a/Project.toml b/Project.toml index ae6968f2..e4364e03 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "ModelPredictiveControl" uuid = "61f9bdb8-6ae4-484a-811f-bbf86720c31c" authors = ["Francis Gagnon"] -version = "1.0.2" +version = "1.1.0" [deps] ControlSystemsBase = "aaaaaaaa-a6ca-5380-bf3e-84a91bcd477e" diff --git a/README.md b/README.md index 4bdcd66c..203b5f8d 100644 --- a/README.md +++ b/README.md @@ -70,71 +70,66 @@ for more detailed examples. ## Features -### Legend - -- [x] implemented feature -- [ ] planned feature - ### Model Predictive Control Features -- [x] linear and nonlinear plant models exploiting multiple dispatch -- [x] model linearization based on automatic differentiation (exact Jacobians) -- [x] supported objective function terms: - - [x] output setpoint tracking - - [x] move suppression - - [x] input setpoint tracking - - [x] terminal costs - - [x] custom economic costs (economic model predictive control) -- [x] adaptive linear model predictive controller - - [x] manual model modification - - [x] automatic successive linearization of a nonlinear model - - [x] objective function weights and covariance matrices modification -- [x] explicit predictive controller for problems without constraint -- [x] online-tunable soft and hard constraints on: - - [x] output predictions - - [x] manipulated inputs - - [x] manipulated inputs increments - - [x] terminal states to ensure nominal stability -- [x] custom economic inequality constraints (soft or hard) -- [x] supported feedback strategy: - - [x] state estimator (see State Estimation features) - - [x] internal model structure with a custom stochastic model -- [x] automatic model augmentation with integrating states for offset-free tracking -- [x] support for unmeasured model outputs -- [x] feedforward action with measured disturbances that supports direct transmission -- [x] custom predictions for: - - [x] output setpoints - - [x] measured disturbances -- [x] easy integration with `Plots.jl` -- [x] optimization based on `JuMP.jl`: - - [x] quickly compare multiple optimizers - - [x] nonlinear solvers relying on automatic differentiation (exact derivative) -- [x] additional information about the optimum to ease troubleshooting -- [x] real-time control loop features: - - [x] implementations that carefully limits the allocations - - [x] simple soft real-time utilities +- linear and nonlinear plant models exploiting multiple dispatch +- model linearization based on automatic differentiation (exact Jacobians) +- supported objective function terms: + - output setpoint tracking + - move suppression + - input setpoint tracking + - terminal costs + - custom economic costs (economic model predictive control) +- adaptive linear model predictive controller + - manual model modification + - automatic successive linearization of a nonlinear model + - objective function weights and covariance matrices modification +- explicit predictive controller for problems without constraint +- online-tunable soft and hard constraints on: + - output predictions + - manipulated inputs + - manipulated inputs increments + - terminal states to ensure nominal stability +- custom nonlinear inequality constraints (soft or hard) +- supported feedback strategy: + - state estimator (see State Estimation features) + - internal model structure with a custom stochastic model +- automatic model augmentation with integrating states for offset-free tracking +- support for unmeasured model outputs +- feedforward action with measured disturbances that supports direct transmission +- custom predictions for: + - output setpoints + - measured disturbances +- easy integration with `Plots.jl` +- optimization based on `JuMP.jl`: + - quickly compare multiple optimizers + - nonlinear solvers relying on automatic differentiation (exact derivative) +- additional information about the optimum to ease troubleshooting +- real-time control loop features: + - implementations that carefully limits the allocations + - simple soft real-time utilities ### State Estimation Features -- [x] supported state estimators/observers: - - [x] steady-state Kalman filter - - [x] Kalman filter - - [x] Luenberger observer - - [x] internal model structure - - [x] extended Kalman filter - - [x] unscented Kalman filter - - [x] moving horizon estimator -- [x] easily estimate unmeasured disturbances by adding one or more integrators at the: - - [x] manipulated inputs - - [x] measured outputs -- [x] bumpless manual to automatic transfer for control with a proper initial estimate -- [x] estimators in two possible forms: - - [x] filter (or current) form to improve accuracy and robustness - - [x] predictor (or delayed) form to reduce computational load -- [x] moving horizon estimator in two formulations: - - [x] linear plant models (quadratic optimization) - - [x] nonlinear plant models (nonlinear optimization) -- [x] moving horizon estimator online-tunable soft and hard constraints on: - - [x] state estimates - - [x] process noise estimates - - [x] sensor noise estimates +- supported state estimators/observers: + - steady-state Kalman filter + - Kalman filter + - Luenberger observer + - internal model structure + - extended Kalman filter + - unscented Kalman filter + - moving horizon estimator +- easily estimate unmeasured disturbances by adding one or more integrators at the: + - manipulated inputs + - measured outputs +- bumpless manual to automatic transfer for control with a proper initial estimate +- estimators in two possible forms: + - filter (or current) form to improve accuracy and robustness + - predictor (or delayed) form to reduce computational load +- moving horizon estimator in two formulations: + - linear plant models (quadratic optimization) + - nonlinear plant models (nonlinear optimization) +- moving horizon estimator online-tunable soft and hard constraints on: + - state estimates + - process noise estimates + - sensor noise estimates diff --git a/docs/src/public/predictive_control.md b/docs/src/public/predictive_control.md index 23cf0501..729d00dd 100644 --- a/docs/src/public/predictive_control.md +++ b/docs/src/public/predictive_control.md @@ -12,8 +12,9 @@ assumes by default that the current model mismatch estimation is constant in the (same approach as dynamic matrix control, DMC). !!! info - The nomenclature uses capital letters for time series (and matrices) and hats for the - predictions (and estimations, for state estimators). + The nomenclature uses boldfaces for vectors or matrices, capital boldface letters for + vectors representing time series (and also for matrices), and hats for the predictions + (and also for the observer estimations). To be precise, at the ``k``th control period, the vectors that encompass the future measured disturbances ``\mathbf{d̂}``, model outputs ``\mathbf{ŷ}`` and setpoints ``\mathbf{r̂_y}`` @@ -31,7 +32,9 @@ over the prediction horizon ``H_p`` are defined as: \end{bmatrix} ``` -The vectors for the manipulated input ``\mathbf{u}`` are shifted by one time step: +in which ``\mathbf{D̂}``, ``\mathbf{Ŷ}`` and ``\mathbf{R̂_y}`` are vectors of `nd*Hp`, `ny*Hp` +and `ny*Hp` elements, respectively. The vectors for the manipulated input ``\mathbf{u}`` +are shifted by one time step: ```math \mathbf{U} = \begin{bmatrix} @@ -42,10 +45,10 @@ The vectors for the manipulated input ``\mathbf{u}`` are shifted by one time ste \end{bmatrix} ``` -Defining the manipulated input increment as ``\mathbf{Δu}(k+j) = -\mathbf{u}(k+j) - \mathbf{u}(k+j-1)``, the control horizon ``H_c`` enforces that -``\mathbf{Δu}(k+j) = \mathbf{0}`` for ``j ≥ H_c``. For this reason, the vector that collects -them is truncated up to ``k+H_c-1``: +in which ``\mathbf{U}`` and ``\mathbf{R̂_u}`` are both vectors of `nu*Hp` elements. Defining +the manipulated input increment as ``\mathbf{Δu}(k+j) = \mathbf{u}(k+j) - \mathbf{u}(k+j-1)``, +the control horizon ``H_c`` enforces that ``\mathbf{Δu}(k+j) = \mathbf{0}`` for ``j ≥ H_c``. +For this reason, the vector that collects them is truncated up to ``k+H_c-1``: ```math \mathbf{ΔU} = @@ -54,6 +57,8 @@ them is truncated up to ``k+H_c-1``: \end{bmatrix} ``` +in which ``\mathbf{ΔU}`` is a vector of `nu*Hc` elements. + ## PredictiveController ```@docs diff --git a/docs/src/public/sim_model.md b/docs/src/public/sim_model.md index 6c13b0eb..a6f09658 100644 --- a/docs/src/public/sim_model.md +++ b/docs/src/public/sim_model.md @@ -11,6 +11,12 @@ simulators by calling [`evaloutput`](@ref) and [`updatestate!`](@ref) methods on the states ``\mathbf{x}`` are stored inside [`SimModel`](@ref) instances. Use [`setstate!`](@ref) method to manually modify them. +!!! info + The nomenclature in this page introduces the model manipulated input ``\mathbf{u}``, + measured disturbances ``\mathbf{d}``, state ``\mathbf{x}`` and output ``\mathbf{y}``, + four vectors of `nu`, `nd`, `nx` and `ny` elements, respectively. The ``\mathbf{z}`` + vector combines the elements of ``\mathbf{u}`` and ``\mathbf{d}``. + ## SimModel ```@docs diff --git a/docs/src/public/state_estim.md b/docs/src/public/state_estim.md index 4a88f56b..2d2d9a79 100644 --- a/docs/src/public/state_estim.md +++ b/docs/src/public/state_estim.md @@ -28,8 +28,10 @@ solving the MPC problem with [`moveinput!`](@ref), for when the estimations are [^1]: also denoted ``\mathbf{x̂}_{k|k}`` [elsewhere](https://en.wikipedia.org/wiki/Kalman_filter). !!! info - All the estimators support measured ``\mathbf{y^m}`` and unmeasured ``\mathbf{y^u}`` - model outputs, where ``\mathbf{y}`` refers to all of them. + The nomenclature in this page introduces the estimated state ``\mathbf{x̂}`` and output + ``\mathbf{ŷ}`` vectors of respectively `nx̂` and `ny` elements. Also, all the estimators + support measured ``\mathbf{y^m}`` (`nym` elements) and unmeasured ``\mathbf{y^u}`` + (`nyu` elements) model output, where ``\mathbf{y}`` refers to all of them. ## StateEstimator diff --git a/src/controller/execute.jl b/src/controller/execute.jl index 69e82fa4..040853da 100644 --- a/src/controller/execute.jl +++ b/src/controller/execute.jl @@ -120,8 +120,8 @@ function getinfo(mpc::PredictiveController{NT}) where NT<:Real Ȳ, Ū = similar(mpc.Yop), similar(mpc.Uop) Ŷe, Ue = Vector{NT}(undef, nŶe), Vector{NT}(undef, nUe) Ŷ0, x̂0end = predict!(Ŷ0, x̂0, x̂0next, u0, û0, mpc, model, mpc.ΔŨ) - Ŷe, Ue = extended_predictions!(Ŷe, Ue, Ū, mpc, model, Ŷ0, mpc.ΔŨ) - J = obj_nonlinprog!(Ȳ, Ū, mpc, model, Ŷe, Ue, mpc.ΔŨ) + Ue, Ŷe = extended_predictions!(Ue, Ŷe, Ū, mpc, model, Ŷ0, mpc.ΔŨ) + J = obj_nonlinprog!(Ȳ, Ū, mpc, model, Ue, Ŷe, mpc.ΔŨ) U = Ū U .= @views Ue[1:end-model.nu] Ŷ = Ȳ @@ -362,28 +362,28 @@ function predict!(Ŷ0, x̂0, x̂0next, u0, û0, mpc::PredictiveController, mod end """ - extended_predictions!(Ŷe, Ue, Ū, mpc, model, Ŷ0, ΔŨ) -> Ŷe, Ue + extended_predictions!(Ue, Ŷe, Ū, mpc, model, Ŷ0, ΔŨ) -> Ŷe, Ue -Compute the extended predictions `Ŷe` and `Ue` for the nonlinear optimization. +Compute the extended vectors `Ue` and `Ŷe` and for the nonlinear optimization. -The function mutates `Ŷe`, `Ue` and `Ū` in arguments, without assuming any initial values. +The function mutates `Ue`, `Ŷe` and `Ū` in arguments, without assuming any initial values. """ -function extended_predictions!(Ŷe, Ue, Ū, mpc, model, Ŷ0, ΔŨ) +function extended_predictions!(Ue, Ŷe, Ū, mpc, model, Ŷ0, ΔŨ) ny, nu = model.ny, model.nu - # --- extended output predictions Ŷe = [ŷ(k); Ŷ] --- - Ŷe[1:ny] .= mpc.ŷ - Ŷe[ny+1:end] .= Ŷ0 .+ mpc.Yop # --- extended manipulated inputs Ue = [U; u(k+Hp-1)] --- U0 = Ū U0 .= mul!(U0, mpc.S̃, ΔŨ) .+ mpc.T_lastu0 Ue[1:end-nu] .= U0 .+ mpc.Uop # u(k + Hp) = u(k + Hp - 1) since Δu(k+Hp) = 0 (because Hc ≤ Hp): Ue[end-nu+1:end] .= @views Ue[end-2nu+1:end-nu] - return Ŷe, Ue + # --- extended output predictions Ŷe = [ŷ(k); Ŷ] --- + Ŷe[1:ny] .= mpc.ŷ + Ŷe[ny+1:end] .= Ŷ0 .+ mpc.Yop + return Ue, Ŷe end """ - obj_nonlinprog!( _ , _ , mpc::PredictiveController, model::LinModel, Ŷe, Ue, ΔŨ) + obj_nonlinprog!( _ , _ , mpc::PredictiveController, model::LinModel, Ue, Ŷe, ΔŨ) Nonlinear programming objective function when `model` is a [`LinModel`](@ref). @@ -392,15 +392,15 @@ also be called on any [`PredictiveController`](@ref)s to evaluate the objective at specific `Ue`, `Ŷe` and `ΔŨ`, values. It does not mutate any argument. """ function obj_nonlinprog!( - _, _, mpc::PredictiveController, model::LinModel, Ŷe, Ue, ΔŨ::AbstractVector{NT} + _, _, mpc::PredictiveController, model::LinModel, Ue, Ŷe, ΔŨ::AbstractVector{NT} ) where NT <: Real JQP = obj_quadprog(ΔŨ, mpc.H̃, mpc.q̃) + mpc.r[] - E_JE = obj_econ!(Ue, Ŷe, mpc, model) + E_JE = obj_econ(mpc, model, Ue, Ŷe) return JQP + E_JE end """ - obj_nonlinprog!(Ȳ, Ū, mpc::PredictiveController, model::SimModel, Ŷe, Ue, ΔŨ) + obj_nonlinprog!(Ȳ, Ū, mpc::PredictiveController, model::SimModel, Ue, Ŷe, ΔŨ) Nonlinear programming objective method when `model` is not a [`LinModel`](@ref). The function `dot(x, A, x)` is a performant way of calculating `x'*A*x`. This method mutates @@ -408,7 +408,7 @@ function `dot(x, A, x)` is a performant way of calculating `x'*A*x`. This method `Ŷe` and `Ue` arguments). """ function obj_nonlinprog!( - Ȳ, Ū, mpc::PredictiveController, model::SimModel, Ŷe, Ue, ΔŨ::AbstractVector{NT} + Ȳ, Ū, mpc::PredictiveController, model::SimModel, Ue, Ŷe, ΔŨ::AbstractVector{NT} ) where NT<:Real nu, ny = model.nu, model.ny # --- output setpoint tracking term --- @@ -426,12 +426,12 @@ function obj_nonlinprog!( JR̂u = 0.0 end # --- economic term --- - E_JE = obj_econ!(Ue, Ŷe, mpc, model) + E_JE = obj_econ(mpc, model, Ue, Ŷe) return JR̂y + JΔŨ + JR̂u + E_JE end "By default, the economic term is zero." -obj_econ!( _ , _ , ::PredictiveController, ::SimModel) = 0.0 +obj_econ(::PredictiveController, ::SimModel, _ , _ ) = 0.0 @doc raw""" optim_objective!(mpc::PredictiveController) -> ΔŨ diff --git a/src/controller/linmpc.jl b/src/controller/linmpc.jl index 770ba49c..b9500a91 100644 --- a/src/controller/linmpc.jl +++ b/src/controller/linmpc.jl @@ -175,6 +175,7 @@ LinMPC controller with a sample time Ts = 4.0 s, OSQP optimizer, SteadyKalmanFil | ``H_p`` | prediction horizon (integer) | `()` | | ``H_c`` | control horizon (integer) | `()` | | ``\mathbf{ΔU}`` | manipulated input increments over ``H_c`` | `(nu*Hc,)` | + | ``\mathbf{D̂}`` | predicted measured disturbances over ``H_p`` | `(nd*Hp,)` | | ``\mathbf{Ŷ}`` | predicted outputs over ``H_p`` | `(ny*Hp,)` | | ``\mathbf{U}`` | manipulated inputs over ``H_p`` | `(nu*Hp,)` | | ``\mathbf{R̂_y}`` | predicted output setpoints over ``H_p`` | `(ny*Hp,)` | diff --git a/src/controller/nonlinmpc.jl b/src/controller/nonlinmpc.jl index 4a0db080..29671e88 100644 --- a/src/controller/nonlinmpc.jl +++ b/src/controller/nonlinmpc.jl @@ -78,6 +78,7 @@ struct NonLinMPC{ # dummy vals (updated just before optimization): d0, D̂0, D̂e = zeros(NT, nd), zeros(NT, nd*Hp), zeros(NT, nd + nd*Hp) Uop, Yop, Dop = repeat(model.uop, Hp), repeat(model.yop, Hp), repeat(model.dop, Hp) + test_custom_functions(NT, model, JE, gc!, nc, Uop, Yop, Dop, p) nΔŨ = size(Ẽ, 2) ΔŨ = zeros(NT, nΔŨ) buffer = PredictiveControllerBuffer{NT}(nu, ny, nd, Hp) @@ -118,30 +119,28 @@ controller minimizes the following objective function at each discrete time ``k` ``` subject to [`setconstraint!`](@ref) bounds, and the custom inequality constraints: ```math -\mathbf{g_c}\Big(\mathbf{U_e}, \mathbf{Ŷ_e}, \mathbf{D̂_e}, \mathbf{p}, ϵ(k)\Big) ≤ \mathbf{0} +\mathbf{g_c}(\mathbf{U_e}, \mathbf{Ŷ_e}, \mathbf{D̂_e}, \mathbf{p}, ϵ) ≤ \mathbf{0} ``` The economic function ``J_E`` can penalizes solutions with high economic costs. Setting all the weights to 0 except ``E`` creates a pure economic model predictive controller (EMPC). As a matter of fact, ``J_E`` can be any nonlinear function to customize the objective, even if there is no economic interpretation to it. The arguments of ``J_E`` and ``\mathbf{g_c}`` include the manipulated inputs, predicted outputs and measured disturbances, extended from -``k`` to ``k + H_p`` (inclusively): +``k`` to ``k + H_p`` (inclusively, see Extended Help for more details): ```math \mathbf{U_e} = \begin{bmatrix} \mathbf{U} \\ \mathbf{u}(k+H_p-1) \end{bmatrix} , \quad \mathbf{Ŷ_e} = \begin{bmatrix} \mathbf{ŷ}(k) \\ \mathbf{Ŷ} \end{bmatrix} , \quad \mathbf{D̂_e} = \begin{bmatrix} \mathbf{d}(k) \\ \mathbf{D̂} \end{bmatrix} ``` -since ``H_c ≤ H_p`` implies that ``\mathbf{Δu}(k+H_p) = \mathbf{0}`` or ``\mathbf{u}(k+H_p)= -\mathbf{u}(k+H_p-1)``. The vector ``\mathbf{D̂}`` comprises the measured disturbance -predictions over ``H_p``. The argument ``\mathbf{p}`` is a custom parameter object of any -type, but use a mutable one if you want to modify it later e.g.: a vector. +The argument ``\mathbf{p}`` is a custom parameter object of any type, but use a mutable one +if you want to modify it later e.g.: a vector. See [`LinMPC`](@ref) Extended Help for the +definition of the other variables. !!! tip Replace any of the arguments of ``J_E`` and ``\mathbf{g_c}`` functions with `_` if not needed (see e.g. the default value of `JE` below). -See [`LinMPC`](@ref) for the definition of the other variables. This method uses the default -state estimator : +This method uses the default state estimator: - if `model` is a [`LinModel`](@ref), a [`SteadyKalmanFilter`](@ref) with default arguments; - else, an [`UnscentedKalmanFilter`](@ref) with default arguments. @@ -197,24 +196,36 @@ NonLinMPC controller with a sample time Ts = 10.0 s, Ipopt optimizer, UnscentedK algebra instead of a `for` loop. This feature can accelerate the optimization, especially for the constraint handling, and is not available in any other package, to my knowledge. + The economic cost ``J_E`` and custom constraint ``\mathbf{g_c}`` functions receive the + extended vectors ``\mathbf{U_e}`` (`nu*Hp+nu` elements), ``\mathbf{Ŷ_e}`` (`ny+ny*Hp` + elements) and ``\mathbf{D̂_e}`` (`nd+nd*Hp` elements) as arguments. They all include the + values from ``k`` to ``k + H_p`` (inclusively). The custom constraint also receives the + slack ``ϵ`` (scalar), which is always zero if `Cwt=Inf`. + + More precisely, the last two time steps in ``\mathbf{U_e}`` are forced to be equal, i.e. + ``\mathbf{u}(k+H_p) = \mathbf{u}(k+H_p-1)``, since ``H_c ≤ H_p`` implies that + ``\mathbf{Δu}(k+H_p) = \mathbf{0}``. The vectors ``\mathbf{ŷ}(k)`` and ``\mathbf{d}(k)`` + are the current state estimator output and measured disturbance, respectively, and + ``\mathbf{Ŷ}`` and ``\mathbf{D̂}``, their respective predictions from ``k+1`` to ``k+H_p``. + If `LHS` represents the result of the left-hand side in the inequality + ``\mathbf{g_c}(\mathbf{U_e}, \mathbf{Ŷ_e}, \mathbf{D̂_e}, \mathbf{p}, ϵ) ≤ \mathbf{0}``, + the function `gc` can be implemented in two possible ways: + + 1. **Non-mutating function** (out-of-place): define it as `gc(Ue, Ŷe, D̂e, p, ϵ) -> LHS`. + This syntax is simple and intuitive but it allocates more memory. + 2. **Mutating function** (in-place): define it as `gc!(LHS, Ue, Ŷe, D̂e, p, ϵ) -> nothing`. + This syntax reduces the allocations and potentially the computational burden as well. + + The keyword argument `nc` is the number of elements in `LHS`, and `gc!`, an alias for + the `gc` argument (both `gc` and `gc!` accepts non-mutating and mutating functions). + The optimization relies on [`JuMP`](https://github.com/jump-dev/JuMP.jl) automatic differentiation (AD) to compute the objective and constraint derivatives. Optimizers generally benefit from exact derivatives like AD. However, the [`NonLinModel`](@ref) state-space functions must be compatible with this feature. See [Automatic differentiation](https://jump.dev/JuMP.jl/stable/manual/nlp/#Automatic-differentiation) for common mistakes when writing these functions. - If `LHS` represents the result of the left-hand side in the inequality - ``\mathbf{g_c}(\mathbf{U_e}, \mathbf{Ŷ_e}, \mathbf{D̂_e}, \mathbf{p}, ϵ) ≤ \mathbf{0}``, - the function `gc` can be implemented in two ways: - - 1. **Non-mutating function** (out-of-place): define it as `gc(Ue, Ŷe, D̂e, p, ϵ) -> LHS`. - This syntax is simple and intuitive but it allocates more memory. - 2. **Mutating function** (in-place): define it as `gc!(LHS, Ue, Ŷe, D̂e, p, ϵ) -> nothing`. - This syntax reduces the allocations and potentially the computational burden as well. - - The keyword argument `nc` is the number of elements in the `LHS` vector, and `gc!`, an - alias for the `gc` argument (both accepts non-mutating and mutating functions). Note - that if `Cwt≠Inf`, the attribute `nlp_scaling_max_gradient` of `Ipopt` is set to + Note that if `Cwt≠Inf`, the attribute `nlp_scaling_max_gradient` of `Ipopt` is set to `10/Cwt` (if not already set), to scale the small values of ``ϵ``. """ function NonLinMPC( @@ -383,14 +394,44 @@ function get_mutating_gc(NT, gc) return gc! end -function test_custom_functions(JE, gc!, uop; Uop, dop, Dop, ΔŨ, p) - # TODO: contunue here (important to guide the user, sim! can be used on NonLinModel, - # but there is no similar function for the custom functions of NonLinMPC) - Ue = [Uop; uop] - D̂e = [dop; Dop] - Ŷ0, x̂0next = - Ŷ0, x̂0end = predict!(Ŷ0, x̂0, x̂0next, u0, û0, mpc, model, mpc.ΔŨ) - JE = JE(Uop, Uop, Dop, p) +""" + test_custom_functions(NT, model::SimModel, JE, gc!, nc, Uop, Yop, Dop, p) + +Test the custom functions `JE` and `gc!` at the operating point `Uop`, `Yop`, `Dop`. + +This function is called at the end of `NonLinMPC` construction. It warns the user if the +custom cost `JE` and constraint `gc!` functions crash at `model` operating points. This +should ease troubleshooting of simple bugs e.g.: the user forgets to set the `nc` argument. +""" +function test_custom_functions(NT, model::SimModel, JE, gc!, nc, Uop, Yop, Dop, p) + uop, dop, yop = model.uop, model.dop, model.yop + Ue, Ŷe, D̂e = [Uop; uop], [yop; Yop], [dop; Dop] + try + val::NT = JE(Ue, Ŷe, D̂e, p) + catch err + @warn( + """ + Calling the JE function with Ue, Ŷe, D̂e arguments fixed at uop=$uop, + yop=$yop, dop=$dop failed with the following stacktrace. Did you forget + to set the keyword argument p? + """, + exception=(err, catch_backtrace()) + ) + end + ϵ, gc = 0, Vector{NT}(undef, nc) + try + gc!(gc, Ue, Ŷe, D̂e, p, ϵ) + catch err + @warn( + """ + Calling the gc function with Ue, Ŷe, D̂e, ϵ arguments fixed at uop=$uop, + yop=$yop, dop=$dop, ϵ=0 failed with the following stacktrace. Did you + forget to set the keyword argument p or nc? + """, + exception=(err, catch_backtrace()) + ) + end + return nothing end """ @@ -438,7 +479,7 @@ function init_optimization!(mpc::NonLinMPC, model::SimModel, optim) @operator(optim, J, nΔŨ, Jfunc) @objective(optim, Min, J(ΔŨvar...)) init_nonlincon!(mpc, model, gfunc) - set_nonlincon!(mpc, model, mpc.optim) #TODO: check if this is really necessary !! + set_nonlincon!(mpc, model, mpc.optim) return nothing end @@ -476,13 +517,13 @@ function get_optim_functions(mpc::NonLinMPC, ::JuMP.GenericModel{JNT}) where JNT Ȳ, Ū = get_tmp(Ȳ_cache, ΔŨ1), get_tmp(Ū_cache, ΔŨ1) x̂0, x̂0next = get_tmp(x̂0_cache, ΔŨ1), get_tmp(x̂0next_cache, ΔŨ1) u0, û0 = get_tmp(u0_cache, ΔŨ1), get_tmp(û0_cache, ΔŨ1) - g, gc = get_tmp(g_cache, ΔŨ1), get_tmp(gc_cache, ΔŨ1) + gc = get_tmp(gc_cache, ΔŨ1) Ŷ0, x̂0end = predict!(Ȳ, x̂0, x̂0next, u0, û0, mpc, model, ΔŨ) - Ŷe, Ue = extended_predictions!(Ŷe, Ue, Ū, mpc, model, Ŷ0, ΔŨ) - ϵ = (nϵ == 1) ? ΔŨ[end] : zero(JNT) # ϵ = 0 if nϵ == 0 (meaning no relaxation) + Ue, Ŷe = extended_predictions!(Ue, Ŷe, Ū, mpc, model, Ŷ0, ΔŨ) + ϵ = (nϵ ≠ 0) ? ΔŨ[end] : zero(T) # ϵ = 0 if nϵ == 0 (meaning no relaxation) mpc.con.gc!(gc, Ue, Ŷe, mpc.D̂e, mpc.p, ϵ) g = con_nonlinprog!(g, mpc, model, x̂0end, Ŷ0, gc, ϵ) - return obj_nonlinprog!(Ȳ, Ū, mpc, model, Ŷe, Ue, ΔŨ)::T + return obj_nonlinprog!(Ȳ, Ū, mpc, model, Ue, Ŷe, ΔŨ)::T end function gfunc_i(i, ΔŨtup::NTuple{N, T}) where {N, T<:Real} ΔŨ1 = ΔŨtup[begin] @@ -495,10 +536,10 @@ function get_optim_functions(mpc::NonLinMPC, ::JuMP.GenericModel{JNT}) where JNT Ȳ, Ū = get_tmp(Ȳ_cache, ΔŨ1), get_tmp(Ū_cache, ΔŨ1) x̂0, x̂0next = get_tmp(x̂0_cache, ΔŨ1), get_tmp(x̂0next_cache, ΔŨ1) u0, û0 = get_tmp(u0_cache, ΔŨ1), get_tmp(û0_cache, ΔŨ1) - g, gc = get_tmp(g_cache, ΔŨ1), get_tmp(gc_cache, ΔŨ1) + gc = get_tmp(gc_cache, ΔŨ1) Ŷ0, x̂0end = predict!(Ȳ, x̂0, x̂0next, u0, û0, mpc, model, ΔŨ) - Ŷe, Ue = extended_predictions!(Ŷe, Ue, Ū, mpc, model, Ŷ0, ΔŨ) - ϵ = (nϵ == 1) ? ΔŨ[end] : zero(JNT) # ϵ = 0 if nϵ == 0 (meaning no relaxation) + Ue, Ŷe = extended_predictions!(Ue, Ŷe, Ū, mpc, model, Ŷ0, ΔŨ) + ϵ = (nϵ ≠ 0) ? ΔŨ[end] : zero(T) # ϵ = 0 if nϵ == 0 (meaning no relaxation) mpc.con.gc!(gc, Ue, Ŷe, mpc.D̂e, mpc.p, ϵ) g = con_nonlinprog!(g, mpc, model, x̂0end, Ŷ0, gc, ϵ) end @@ -656,7 +697,7 @@ function con_nonlinprog!(g, mpc::NonLinMPC, ::SimModel, x̂0end, Ŷ0, gc, ϵ) end "Evaluate the economic term `E*JE` of the objective function for [`NonLinMPC`](@ref)." -function obj_econ!(Ue, Ŷe, mpc::NonLinMPC, model::SimModel) +function obj_econ(mpc::NonLinMPC, model::SimModel, Ue, Ŷe) E_JE = iszero(mpc.E) ? 0.0 : mpc.E*mpc.JE(Ue, Ŷe, mpc.D̂e, mpc.p) return E_JE end \ No newline at end of file diff --git a/src/estimator/mhe/construct.jl b/src/estimator/mhe/construct.jl index ca85cb8d..6c092250 100644 --- a/src/estimator/mhe/construct.jl +++ b/src/estimator/mhe/construct.jl @@ -1298,7 +1298,7 @@ function get_optim_functions( estim::MovingHorizonEstimator, ::JuMP.GenericModel{JNT} ) where {JNT <: Real} model, con = estim.model, estim.con - nx̂, nym, nŷ, nu, He = estim.nx̂, estim.nym, model.ny, model.nu, estim.He + nx̂, nym, nŷ, nu, nϵ, He = estim.nx̂, estim.nym, model.ny, model.nu, estim.nϵ, estim.He nV̂, nX̂, ng, nZ̃ = He*nym, He*nx̂, length(con.i_g), length(estim.Z̃) Nc = nZ̃ + 3 Z̃_cache::DiffCache{Vector{JNT}, Vector{JNT}} = DiffCache(zeros(JNT, nZ̃), Nc) @@ -1317,8 +1317,8 @@ function get_optim_functions( V̂, X̂0 = get_tmp(V̂_cache, Z̃1), get_tmp(X̂0_cache, Z̃1) û0, ŷ0 = get_tmp(û0_cache, Z̃1), get_tmp(ŷ0_cache, Z̃1) V̂, X̂0 = predict!(V̂, X̂0, û0, ŷ0, estim, model, Z̃) - g = get_tmp(g_cache, Z̃1) - g = con_nonlinprog!(g, estim, model, X̂0, V̂, Z̃) + ϵ = (nϵ ≠ 0) ? Z̃[begin] : zero(T) # ϵ = 0 if Cwt=Inf (meaning: no relaxation) + g = con_nonlinprog!(g, estim, model, X̂0, V̂, ϵ) x̄ = get_tmp(x̄_cache, Z̃1) return obj_nonlinprog!(x̄, estim, model, V̂, Z̃)::T end @@ -1332,7 +1332,8 @@ function get_optim_functions( Z̃[i] = Z̃tup[i] # Z̃ .= Z̃tup seems to produce a type instability end V̂, X̂0 = predict!(V̂, X̂0, û0, ŷ0, estim, model, Z̃) - g = con_nonlinprog!(g, estim, model, X̂0, V̂, Z̃) + ϵ = (nϵ ≠ 0) ? Z̃[begin] : zero(T) # ϵ = 0 if Cwt=Inf (meaning: no relaxation) + g = con_nonlinprog!(g, estim, model, X̂0, V̂, ϵ) end return g[i] end diff --git a/src/estimator/mhe/execute.jl b/src/estimator/mhe/execute.jl index 3e29abc1..29a5ec00 100644 --- a/src/estimator/mhe/execute.jl +++ b/src/estimator/mhe/execute.jl @@ -541,14 +541,13 @@ function predict!(V̂, X̂0, û0, ŷ0, estim::MovingHorizonEstimator, model::S end """ - con_nonlinprog!(g, estim::MovingHorizonEstimator, model::SimModel, X̂0, V̂, Z̃) + con_nonlinprog!(g, estim::MovingHorizonEstimator, model::SimModel, X̂0, V̂, ϵ) Nonlinear constrains for [`MovingHorizonEstimator`](@ref). """ -function con_nonlinprog!(g, estim::MovingHorizonEstimator, ::SimModel, X̂0, V̂, Z̃) +function con_nonlinprog!(g, estim::MovingHorizonEstimator, ::SimModel, X̂0, V̂, ϵ) nX̂con, nX̂ = length(estim.con.X̂0min), estim.nx̂ *estim.Nk[] nV̂con, nV̂ = length(estim.con.V̂min), estim.nym*estim.Nk[] - ϵ = estim.nϵ ≠ 0 ? Z̃[begin] : 0 # ϵ = 0 if Cwt=Inf (meaning: no relaxation) for i in eachindex(g) estim.con.i_g[i] || continue if i ≤ nX̂con diff --git a/test/test_predictive_control.jl b/test/test_predictive_control.jl index 6cbb20f8..eb6579cb 100644 --- a/test/test_predictive_control.jl +++ b/test/test_predictive_control.jl @@ -18,7 +18,10 @@ sys = [ tf(1.90,[1800.0,1]) tf(1.90,[1800.0,1]) tf(1.90,[1800.0,1]); @test mpc5.Ñ_Hc ≈ Diagonal(diagm([repeat(Float64[3, 4], 5); [1e3]])) mpc6 = LinMPC(model, Lwt=[0,1], Hp=15) @test mpc6.L_Hp ≈ Diagonal(diagm(repeat(Float64[0, 1], 15))) - mpc7 = LinMPC(model, optim=JuMP.Model(DAQP.Optimizer)) + mpc7 = @test_logs( + (:warn, "Solving time limit is not supported by the optimizer."), + LinMPC(model, optim=JuMP.Model(DAQP.Optimizer)) + ) @test solver_name(mpc7.optim) == "DAQP" kf = KalmanFilter(model) mpc8 = LinMPC(kf) @@ -37,7 +40,12 @@ sys = [ tf(1.90,[1800.0,1]) tf(1.90,[1800.0,1]) tf(1.90,[1800.0,1]); @test isa(mpc13, LinMPC{Float32}) @test isa(mpc13.optim, JuMP.GenericModel{Float64}) # OSQP does not support Float32 - @test_throws ArgumentError LinMPC(model, Hp=0) + @test_logs( + (:warn, + "prediction horizon Hp (0) ≤ estimated number of delays in model (0), the "* + "closed-loop system may be unstable or zero-gain (unresponsive)"), + @test_throws ArgumentError LinMPC(model, Hp=0) + ) @test_throws ArgumentError LinMPC(model, Hc=0) @test_throws ArgumentError LinMPC(model, Hp=1, Hc=2) @test_throws ArgumentError LinMPC(model, Mwt=[1]) @@ -134,8 +142,11 @@ end preparestate!(mpc1, [50, 30]) updatestate!(mpc1, mpc1.estim.model.uop, [50, 30]) @test mpc1.estim.x̂0 ≈ [0,0,0,0] - # do not call preparestate! before moveinput! for the warning: - moveinput!(mpc1, [10, 50]) + @test_logs( + (:warn, "preparestate! should be called before moveinput! with current estimators"), + (:warn, "preparestate! should be called before evaloutput with current estimators"), + moveinput!(mpc1, [10, 50]) + ) @test_throws ArgumentError updatestate!(mpc1, [0,0]) end @@ -219,6 +230,10 @@ end info = getinfo(mpc) @test info[:ΔU][begin] ≈ -1.5 atol=1e-1 @test info[:U][end] ≈ -3 atol=1e-1 + moveinput!(mpc, [10]) + info = getinfo(mpc) + @test info[:ΔU][begin] ≈ 1.5 atol=1e-1 + @test info[:U][end] ≈ 3 atol=1e-1 setconstraint!(mpc, umin=[-10], umax=[10]) setconstraint!(mpc, Δumin=[-15], Δumax=[15]) @@ -226,6 +241,9 @@ end moveinput!(mpc, [-10]) info = getinfo(mpc) @test info[:Ŷ][end] ≈ -0.5 atol=1e-1 + moveinput!(mpc, [10]) + info = getinfo(mpc) + @test info[:Ŷ][end] ≈ 0.5 atol=1e-1 setconstraint!(mpc, umin=[-10], umax=[10]) setconstraint!(mpc, Δumin=[-15], Δumax=[15]) @@ -234,6 +252,10 @@ end info = getinfo(mpc) @test info[:Ŷ][end] ≈ -10 atol=1e-1 @test info[:Ŷ][begin] ≈ -0.5 atol=1e-1 + moveinput!(mpc, [10]) + info = getinfo(mpc) + @test info[:Ŷ][end] ≈ 10 atol=1e-1 + @test info[:Ŷ][begin] ≈ 0.5 atol=1e-1 setconstraint!(mpc, umin=[-1e3], umax=[+1e3]) setconstraint!(mpc, Δumin=[-1e3], Δumax=[+1e3]) @@ -242,6 +264,9 @@ end moveinput!(mpc, [-10]) info = getinfo(mpc) @test info[:x̂end][1] ≈ 0 atol=1e-1 + moveinput!(mpc, [10]) + info = getinfo(mpc) + @test info[:x̂end][1] ≈ 0 atol=1e-1 end @testset "LinMPC terminal cost" begin @@ -491,9 +516,9 @@ end @test nmpc5.Ñ_Hc ≈ Diagonal(diagm([repeat(Float64[3, 4], 5); [1e3]])) nmpc6 = NonLinMPC(nonlinmodel, Hp=15, Lwt=[0,1]) @test nmpc6.L_Hp ≈ Diagonal(diagm(repeat(Float64[0, 1], 15))) - nmpc7 = NonLinMPC(nonlinmodel, Hp=15, Ewt=1e-3, JE=(Ue,Ŷe,D̂e,p) -> p*Ue.*Ŷe.*D̂e, p=2) + nmpc7 = NonLinMPC(nonlinmodel, Hp=15, Ewt=1e-3, JE=(Ue,Ŷe,D̂e,p) -> p*dot(Ue,Ŷe)+sum(D̂e), p=10) @test nmpc7.E == 1e-3 - @test nmpc7.JE([1,2],[3,4],[4,6],2) == 2*[1,2].*[3,4].*[4,6] + @test nmpc7.JE([1,2],[3,4],[4,6],10) == 10*dot([1,2],[3,4])+sum([4,6]) optim = JuMP.Model(optimizer_with_attributes(Ipopt.Optimizer, "nlp_scaling_max_gradient"=>1.0)) nmpc8 = NonLinMPC(nonlinmodel, Hp=15, optim=optim) @test solver_name(nmpc8.optim) == "Ipopt" @@ -513,6 +538,10 @@ end @test nmpc13.Ñ_Hc ≈ Diagonal([0.1,0.11,0.12,0.13]) nmcp14 = NonLinMPC(nonlinmodel, Hp=10, L_Hp=Diagonal(collect(0.001:0.001:0.02))) @test nmcp14.L_Hp ≈ Diagonal(collect(0.001:0.001:0.02)) + nmpc15 = NonLinMPC(nonlinmodel, Hp=10, gc=(Ue,Ŷe,D̂e,p,ϵ)-> [p*dot(Ue,Ŷe)+sum(D̂e)+ϵ], nc=1, p=10) + LHS = zeros(1) + nmpc15.con.gc!(LHS,[1,2],[3,4],[4,6],10,0.1) + @test LHS ≈ [10*dot([1,2],[3,4])+sum([4,6])+0.1] nonlinmodel2 = NonLinModel{Float32}(f, h, Ts, 2, 4, 2, 1, solver=nothing) nmpc15 = NonLinMPC(nonlinmodel2, Hp=15) @@ -521,7 +550,12 @@ end @test_throws ArgumentError NonLinMPC(nonlinmodel, Hp=15, Ewt=[1, 1]) @test_throws ArgumentError NonLinMPC(nonlinmodel) - @test_throws ErrorException NonLinMPC(nonlinmodel, Hp=15, JE=(_,_,_)->0.0) + @test_throws ErrorException NonLinMPC(nonlinmodel, Hp=15, JE = (_,_,_)->0.0) + @test_throws ErrorException NonLinMPC(nonlinmodel, Hp=15, gc = (_,_,_,_)->[0.0], nc=1) + @test_throws ErrorException NonLinMPC(nonlinmodel, Hp=15, gc! = (_,_,_,_)->[0.0], nc=1) + + @test_logs (:warn, Regex(".*")) NonLinMPC(nonlinmodel, Hp=15, JE=(Ue,_,_,_)->Ue) + @test_logs (:warn, Regex(".*")) NonLinMPC(nonlinmodel, Hp=15, gc=(Ue,_,_,_,_)->Ue, nc=0) end @testset "NonLinMPC moves and getinfo" begin @@ -680,8 +714,10 @@ end end @testset "NonLinMPC constraint violation" begin + gc( _ , Ŷe, _ , p , ϵ) = p[]*(Ŷe .- 3.14 .- ϵ) + linmodel = LinModel(tf([2], [10000, 1]), 3000.0) - nmpc_lin = NonLinMPC(linmodel, Hp=50, Hc=5) + nmpc_lin = NonLinMPC(linmodel, Hp=50, Hc=5, gc=gc, nc=50+1, p=[0]) setconstraint!(nmpc_lin, x̂min=[-1e3,-Inf], x̂max=[1e3,+Inf]) setconstraint!(nmpc_lin, umin=[-3], umax=[3]) @@ -690,23 +726,34 @@ end preparestate!(nmpc_lin, [0]) moveinput!(nmpc_lin, [-20]) info = getinfo(nmpc_lin) - @test info[:ΔU][begin] ≈ -1.5 atol=1e-2 - @test info[:U][end] ≈ -3 atol=1e-2 + @test info[:ΔU][begin] ≈ -1.5 atol=1e-1 + @test info[:U][end] ≈ -3 atol=1e-1 + moveinput!(nmpc_lin, [20]) + info = getinfo(nmpc_lin) + @test info[:ΔU][begin] ≈ 1.5 atol=1e-1 + @test info[:U][end] ≈ 3 atol=1e-1 setconstraint!(nmpc_lin, umin=[-10], umax=[10]) setconstraint!(nmpc_lin, Δumin=[-15], Δumax=[15]) setconstraint!(nmpc_lin, ymin=[-0.5], ymax=[0.5]) moveinput!(nmpc_lin, [-20]) info = getinfo(nmpc_lin) - @test info[:Ŷ][end] ≈ -0.5 atol=1e-2 + @test info[:Ŷ][end] ≈ -0.5 atol=1e-1 + moveinput!(nmpc_lin, [20]) + info = getinfo(nmpc_lin) + @test info[:Ŷ][end] ≈ 0.5 atol=1e-1 setconstraint!(nmpc_lin, umin=[-10], umax=[10]) setconstraint!(nmpc_lin, Δumin=[-15], Δumax=[15]) setconstraint!(nmpc_lin, Ymin=[-0.5; fill(-100, 49)], Ymax=[0.5; fill(+100, 49)]) moveinput!(nmpc_lin, [-10]) info = getinfo(nmpc_lin) - @test info[:Ŷ][end] ≈ -10 atol=1e-2 - @test info[:Ŷ][begin] ≈ -0.5 atol=1e-2 + @test info[:Ŷ][end] ≈ -10 atol=1e-1 + @test info[:Ŷ][begin] ≈ -0.5 atol=1e-1 + moveinput!(nmpc_lin, [10]) + info = getinfo(nmpc_lin) + @test info[:Ŷ][end] ≈ 10 atol=1e-1 + @test info[:Ŷ][begin] ≈ 0.5 atol=1e-1 setconstraint!(nmpc_lin, umin=[-1e3], umax=[+1e3]) setconstraint!(nmpc_lin, Δumin=[-1e3], Δumax=[+1e3]) @@ -715,11 +762,24 @@ end moveinput!(nmpc_lin, [-10]) info = getinfo(nmpc_lin) @test info[:x̂end][1] ≈ 0 atol=1e-1 + moveinput!(nmpc_lin, [10]) + info = getinfo(nmpc_lin) + @test info[:x̂end][1] ≈ 0 atol=1e-1 + + nmpc_lin.p[] = 1 + setconstraint!(nmpc_lin, x̂min=[-1e3,-Inf], x̂max=[1e3,+Inf]) + setconstraint!(nmpc_lin, umin=[-10], umax=[10]) + setconstraint!(nmpc_lin, Δumin=[-15], Δumax=[15]) + setconstraint!(nmpc_lin, ymin=[-100], ymax=[100]) + moveinput!(nmpc_lin, [20]) + info = getinfo(nmpc_lin) + @test info[:Ŷ][end] ≈ 3.14 atol=1e-1 + @test info[:Ŷ][begin] ≈ 3.14 atol=1e-1 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) + nmpc = NonLinMPC(nonlinmodel, Hp=50, Hc=5, gc=gc, nc=50+1, p=[0]) setconstraint!(nmpc, x̂min=[-1e3,-Inf], x̂max=[1e3,+Inf]) setconstraint!(nmpc, umin=[-3], umax=[3]) @@ -728,23 +788,34 @@ end preparestate!(nmpc, [0]) moveinput!(nmpc, [-20]) info = getinfo(nmpc) - @test info[:ΔU][begin] ≈ -1.5 atol=1e-2 - @test info[:U][end] ≈ -3 atol=1e-2 + @test info[:ΔU][begin] ≈ -1.5 atol=1e-1 + @test info[:U][end] ≈ -3 atol=1e-1 + moveinput!(nmpc, [20]) + info = getinfo(nmpc) + @test info[:ΔU][begin] ≈ 1.5 atol=1e-1 + @test info[:U][end] ≈ 3 atol=1e-1 setconstraint!(nmpc, umin=[-10], umax=[10]) setconstraint!(nmpc, Δumin=[-15], Δumax=[15]) setconstraint!(nmpc, ymin=[-0.5], ymax=[0.5]) moveinput!(nmpc, [-20]) info = getinfo(nmpc) - @test info[:Ŷ][end] ≈ -0.5 atol=1e-2 + @test info[:Ŷ][end] ≈ -0.5 atol=1e-1 + moveinput!(nmpc, [20]) + info = getinfo(nmpc) + @test info[:Ŷ][end] ≈ 0.5 atol=1e-1 setconstraint!(nmpc, umin=[-10], umax=[10]) setconstraint!(nmpc, Δumin=[-15], Δumax=[15]) setconstraint!(nmpc, Ymin=[-0.5; fill(-100, 49)], Ymax=[0.5; fill(+100, 49)]) moveinput!(nmpc, [-10]) info = getinfo(nmpc) - @test info[:Ŷ][end] ≈ -10 atol=1e-2 - @test info[:Ŷ][begin] ≈ -0.5 atol=1e-2 + @test info[:Ŷ][end] ≈ -10 atol=1e-1 + @test info[:Ŷ][begin] ≈ -0.5 atol=1e-1 + moveinput!(nmpc, [10]) + info = getinfo(nmpc) + @test info[:Ŷ][end] ≈ 10 atol=1e-1 + @test info[:Ŷ][begin] ≈ 0.5 atol=1e-1 setconstraint!(nmpc, umin=[-1e3], umax=[+1e3]) setconstraint!(nmpc, Δumin=[-1e3], Δumax=[+1e3]) @@ -753,6 +824,19 @@ end moveinput!(nmpc, [-10]) info = getinfo(nmpc) @test info[:x̂end][1] ≈ 0 atol=1e-1 + moveinput!(nmpc, [10]) + info = getinfo(nmpc) + @test info[:x̂end][1] ≈ 0 atol=1e-1 + + nmpc.p[] = 1 + setconstraint!(nmpc, x̂min=[-1e3,-Inf], x̂max=[1e3,+Inf]) + setconstraint!(nmpc, umin=[-10], umax=[10]) + setconstraint!(nmpc, Δumin=[-15], Δumax=[15]) + setconstraint!(nmpc, ymin=[-100], ymax=[100]) + moveinput!(nmpc, [20]) + info = getinfo(nmpc) + @test info[:Ŷ][end] ≈ 3.14 atol=1e-1 + @test info[:Ŷ][begin] ≈ 3.14 atol=1e-1 end @testset "NonLinMPC set model" begin