Merge pull request #50 from JuliaControl/clean_namespace

cleaning the namespace with `import`s instead of `using`s

Author: franckgaga

src/ModelPredictiveControl.jl

@@ -1,15 +1,25 @@
 module ModelPredictiveControl
 
 using PrecompileTools
+using LinearAlgebra
+using Random: randn
 
-using LinearAlgebra, Random
 using RecipesBase
-using ControlSystemsBase
 using ForwardDiff
-using JuMP
-using PreallocationTools
-import OSQP, Ipopt
 
+import ControlSystemsBase
+import ControlSystemsBase: ss, tf, delay
+import ControlSystemsBase: Continuous, Discrete
+import ControlSystemsBase: StateSpace, TransferFunction, DelayLtiSystem, LTISystem
+import ControlSystemsBase: iscontinuous, isdiscrete, sminreal, minreal, c2d, d2c
+
+import JuMP
+import JuMP: MOIU, MOI, GenericModel, Model, optimizer_with_attributes, register
+import JuMP: @variable, @constraint, @objective, @NLconstraint, @NLobjective
+
+import PreallocationTools: DiffCache, get_tmp
+
+import OSQP, Ipopt
 
 export SimModel, LinModel, NonLinModel
 export DiffSolver, RungeKutta
@@ -28,8 +38,8 @@ include("predictive_control.jl")
 include("plot_sim.jl")
 
 @setup_workload begin
-    # Putting some things in `@setup_workload` instead of `@compile_workload` can reduce the size of the
-    # precompile file and potentially make loading faster.
+    # Putting some things in `@setup_workload` instead of `@compile_workload` can reduce the
+    # size of the precompile file and potentially make loading faster.
     @compile_workload begin
         # all calls in this block will be precompiled, regardless of whether
         # they belong to your package or not (on Julia 1.8 and higher)

src/controller/construct.jl

@@ -130,8 +130,8 @@ function setconstraint!(
 )
     model, con, optim = mpc.estim.model, mpc.con, mpc.optim
     nu, ny, nx̂, Hp, Hc = model.nu, model.ny, mpc.estim.nx̂, mpc.Hp, mpc.Hc
-    notSolvedYet = (termination_status(optim) == OPTIMIZE_NOT_CALLED)
-    C, E, ex̂ = mpc.C, mpc.Ẽ[:, 1:nu*Hc], mpc.con.ẽx̂[:, 1:nu*Hc]
+    notSolvedYet = (JuMP.termination_status(optim) == JuMP.OPTIMIZE_NOT_CALLED)
+    C = mpc.C
     isnothing(Umin)     && !isnothing(umin)     && (Umin    = repeat(umin,    Hp))
     isnothing(Umax)     && !isnothing(umax)     && (Umax    = repeat(umax,    Hp))
     isnothing(ΔUmin)    && !isnothing(Δumin)    && (ΔUmin   = repeat(Δumin,   Hc))
@@ -238,8 +238,8 @@ function setconstraint!(
         A = con.A[con.i_b, :]
         b = con.b[con.i_b]
         ΔŨvar = optim[:ΔŨvar]
-        delete(optim, optim[:linconstraint])
-        unregister(optim, :linconstraint)
+        JuMP.delete(optim, optim[:linconstraint])
+        JuMP.unregister(optim, :linconstraint)
         @constraint(optim, linconstraint, A*ΔŨvar .≤ b)
         setnonlincon!(mpc, model)
     else

src/controller/execute.jl

@@ -134,7 +134,7 @@ end
 By default, add the solution summary `:sol` that can be printed to `info`.
 """
 function addinfo!(info, mpc::PredictiveController)
-    info[:sol] = solution_summary(mpc.optim, verbose=true)
+    info[:sol] = JuMP.solution_summary(mpc.optim, verbose=true)
     return info
 end
 
@@ -280,7 +280,7 @@ function linconstraint!(mpc::PredictiveController, model::LinModel)
     mpc.con.b[(n+1):(n+nx̂)]  .= @. +mpc.con.x̂max - fx̂
     if any(mpc.con.i_b) 
         lincon = mpc.optim[:linconstraint]
-        set_normalized_rhs(lincon, mpc.con.b[mpc.con.i_b])
+        JuMP.set_normalized_rhs(lincon, mpc.con.b[mpc.con.i_b])
     end
     return nothing
 end
@@ -298,7 +298,7 @@ function linconstraint!(mpc::PredictiveController, ::SimModel)
     mpc.con.b[(n+1):(n+nΔŨ)] .= @. +mpc.con.ΔŨmax
     if any(mpc.con.i_b) 
         lincon = mpc.optim[:linconstraint]
-        set_normalized_rhs(lincon, mpc.con.b[mpc.con.i_b])
+        JuMP.set_normalized_rhs(lincon, mpc.con.b[mpc.con.i_b])
     end
     return nothing
 end
@@ -428,38 +428,38 @@ warm-start value.
 function optim_objective!(mpc::PredictiveController{NT}) where {NT<:Real}
     optim = mpc.optim
     model = mpc.estim.model
-    ΔŨvar::Vector{VariableRef} = optim[:ΔŨvar]
+    ΔŨvar::Vector{JuMP.VariableRef} = optim[:ΔŨvar]
     # initial ΔŨ (warm-start): [Δu_{k-1}(k); Δu_{k-1}(k+1); ... ; 0_{nu × 1}; ϵ_{k-1}]
     ϵ0  = !isinf(mpc.C) ? mpc.ΔŨ[end] : empty(mpc.ΔŨ)
     ΔŨ0 = [mpc.ΔŨ[(model.nu+1):(mpc.Hc*model.nu)]; zeros(NT, model.nu); ϵ0]
-    set_start_value.(ΔŨvar, ΔŨ0)
+    JuMP.set_start_value.(ΔŨvar, ΔŨ0)
     set_objective_linear_coef!(mpc, ΔŨvar)
     try
-        optimize!(optim)
+        JuMP.optimize!(optim)
     catch err
         if isa(err, MOI.UnsupportedAttribute{MOI.VariablePrimalStart})
             # reset_optimizer to unset warm-start, set_start_value.(nothing) seems buggy
             MOIU.reset_optimizer(optim)
-            optimize!(optim)
+            JuMP.optimize!(optim)
         else
             rethrow(err)
         end
     end
     if !issolved(optim)
-        status = termination_status(optim)
+        status = JuMP.termination_status(optim)
         if iserror(optim)
             @error("MPC terminated without solution: returning last solution shifted", 
                    status)
         else
             @warn("MPC termination status not OPTIMAL or LOCALLY_SOLVED: keeping "*
                   "solution anyway", status)
         end
-        @debug solution_summary(optim, verbose=true)
+        @debug JuMP.solution_summary(optim, verbose=true)
     end
     if iserror(optim)
         mpc.ΔŨ .= ΔŨ0
     else
-        mpc.ΔŨ .= value.(ΔŨvar)
+        mpc.ΔŨ .= JuMP.value.(ΔŨvar)
     end
     return mpc.ΔŨ
 end

src/controller/linmpc.jl

@@ -243,7 +243,7 @@ function init_optimization!(mpc::LinMPC, optim::JuMP.GenericModel)
     # --- variables and linear constraints ---
     con = mpc.con
     nΔŨ = length(mpc.ΔŨ)
-    set_silent(optim)
+    JuMP.set_silent(optim)
     limit_solve_time(mpc.optim, mpc.estim.model.Ts)
     @variable(optim, ΔŨvar[1:nΔŨ])
     A = con.A[con.i_b, :]
@@ -256,6 +256,6 @@ end
 
 "For [`LinMPC`](@ref), set the QP linear coefficient `q̃` just before optimization."
 function set_objective_linear_coef!(mpc::LinMPC, ΔŨvar)
-    set_objective_coefficient(mpc.optim, ΔŨvar, mpc.q̃)
+    JuMP.set_objective_coefficient(mpc.optim, ΔŨvar, mpc.q̃)
     return nothing
 end

src/controller/nonlinmpc.jl

@@ -272,7 +272,7 @@ function addinfo!(info, mpc::NonLinMPC)
     ŶE = [ŷ; Ŷ]
     D̂E = [d; D̂]
     info[:JE]  = mpc.JE(UE, ŶE, D̂E)
-    info[:sol] = solution_summary(mpc.optim, verbose=true)
+    info[:sol] = JuMP.solution_summary(mpc.optim, verbose=true)
     return info
 end
 
@@ -285,19 +285,19 @@ function init_optimization!(mpc::NonLinMPC, optim::JuMP.GenericModel{JNT}) where
     # --- variables and linear constraints ---
     C, con = mpc.C, mpc.con
     nΔŨ = length(mpc.ΔŨ)
-    set_silent(optim)
+    JuMP.set_silent(optim)
     limit_solve_time(mpc.optim, mpc.estim.model.Ts)
     @variable(optim, ΔŨvar[1:nΔŨ])
     A = con.A[con.i_b, :]
     b = con.b[con.i_b]
     @constraint(optim, linconstraint, A*ΔŨvar .≤ b)
     # --- nonlinear optimization init ---
-    if !isinf(C) && solver_name(optim) == "Ipopt"
+    if !isinf(C) && JuMP.solver_name(optim) == "Ipopt"
         try
-            get_attribute(optim, "nlp_scaling_max_gradient")
+            JuMP.get_attribute(optim, "nlp_scaling_max_gradient")
         catch
             # default "nlp_scaling_max_gradient" to `10.0/C` if not already set:
-            set_attribute(optim, "nlp_scaling_max_gradient", 10.0/C)
+            JuMP.set_attribute(optim, "nlp_scaling_max_gradient", 10.0/C)
         end
     end
     model = mpc.estim.model
@@ -387,22 +387,22 @@ function setnonlincon!(mpc::NonLinMPC, ::NonLinModel)
     optim = mpc.optim
     ΔŨvar = optim[:ΔŨvar]
     con = mpc.con
-    map(con -> delete(optim, con), all_nonlinear_constraints(optim))
+    map(con -> JuMP.delete(optim, con), JuMP.all_nonlinear_constraints(optim))
     for i in findall(.!isinf.(con.Ymin))
         f_sym = Symbol("g_Ymin_$(i)")
-        add_nonlinear_constraint(optim, :($(f_sym)($(ΔŨvar...)) <= 0))
+        JuMP.add_nonlinear_constraint(optim, :($(f_sym)($(ΔŨvar...)) <= 0))
     end
     for i in findall(.!isinf.(con.Ymax))
         f_sym = Symbol("g_Ymax_$(i)")
-        add_nonlinear_constraint(optim, :($(f_sym)($(ΔŨvar...)) <= 0))
+        JuMP.add_nonlinear_constraint(optim, :($(f_sym)($(ΔŨvar...)) <= 0))
     end
     for i in findall(.!isinf.(con.x̂min))
         f_sym = Symbol("g_x̂min_$(i)")
-        add_nonlinear_constraint(optim, :($(f_sym)($(ΔŨvar...)) <= 0))
+        JuMP.add_nonlinear_constraint(optim, :($(f_sym)($(ΔŨvar...)) <= 0))
     end
     for i in findall(.!isinf.(con.x̂max))
         f_sym = Symbol("g_x̂max_$(i)")
-        add_nonlinear_constraint(optim, :($(f_sym)($(ΔŨvar...)) <= 0))
+        JuMP.add_nonlinear_constraint(optim, :($(f_sym)($(ΔŨvar...)) <= 0))
     end
     return nothing
 end

src/estimator/construct.jl

@@ -126,7 +126,7 @@ function augment_model(model::LinModel{NT}, As, Cs_u, Cs_y; verify_obsv=true) wh
     B̂d  = [model.Bd; zeros(NT, nxs, nd)]
     D̂d  = model.Dd
     # observability on Ĉ instead of Ĉm, since it would always return false when nym ≠ ny:
-    if verify_obsv && !observability(Â, Ĉ)[:isobservable]
+    if verify_obsv && !ControlSystemsBase.observability(Â, Ĉ)[:isobservable]
         error("The augmented model is unobservable. You may try to use 0 integrator on "*
               "model integrating outputs with nint_ym parameter. Adding integrators at both "*
               "inputs (nint_u) and outputs (nint_ym) can also violate observability.")
@@ -175,7 +175,7 @@ function default_nint(model::LinModel, i_ym=1:model.ny, nint_u=0)
         As, Cs_u, Cs_y = init_estimstoch(model, i_ym, nint_u, nint_ym)
         Â, _ , Ĉ = augment_model(model, As, Cs_u, Cs_y, verify_obsv=false)
         # observability on Ĉ instead of Ĉm, since it would always return false when nym ≠ ny
-        observability(Â, Ĉ)[:isobservable] || (nint_ym[i] = 0)
+        ControlSystemsBase.observability(Â, Ĉ)[:isobservable] || (nint_ym[i] = 0)
     end
     return nint_ym
 end

src/estimator/kalman.jl

@@ -35,7 +35,7 @@ struct SteadyKalmanFilter{NT<:Real, SM<:LinModel} <: StateEstimator{NT}
             Q̂_kalman = Matrix(Q̂) # Matrix() required for Julia 1.6
             R̂_kalman = zeros(NT, model.ny, model.ny)
             R̂_kalman[i_ym, i_ym] = R̂
-            kalman(Discrete, Â, Ĉ, Q̂_kalman, R̂_kalman)[:, i_ym] 
+            ControlSystemsBase.kalman(Discrete, Â, Ĉ, Q̂_kalman, R̂_kalman)[:, i_ym] 
         catch my_error
             if isa(my_error, ErrorException)
                 error("Cannot compute the optimal Kalman gain K for the "* 

src/estimator/luenberger.jl

@@ -30,7 +30,7 @@ struct Luenberger{NT<:Real, SM<:LinModel} <: StateEstimator{NT}
         nx̂  = model.nx + nxs
         Â, B̂u, Ĉ, B̂d, D̂d = augment_model(model, As, Cs_u, Cs_y)
         K̂ = try
-            place(Â, Ĉ, p̂, :o)[:, i_ym]
+            ControlSystemsBase.place(Â, Ĉ, p̂, :o)[:, i_ym]
         catch
             error("Cannot compute the Luenberger gain K̂ with specified poles p̂.")
         end

src/estimator/mhe.jl

@@ -6,7 +6,7 @@ function Base.show(io::IO, estim::MovingHorizonEstimator)
     nx̂, nym, nyu = estim.nx̂, estim.nym, estim.nyu
     n = maximum(ndigits.((nu, nx̂, nym, nyu, nd))) + 1
     println(io, "$(typeof(estim).name.name) estimator with a sample time "*
-                "Ts = $(estim.model.Ts) s, $(solver_name(estim.optim)) optimizer, "*
+                "Ts = $(estim.model.Ts) s, $(JuMP.solver_name(estim.optim)) optimizer, "*
                 "$(typeof(estim.model).name.name) and:")
     print_estim_dim(io, estim, n)
 end

src/estimator/mhe/construct.jl

@@ -412,7 +412,7 @@ function setconstraint!(
     model, optim, con = estim.model, estim.optim, estim.con
     nx̂, nŵ, nym, He = estim.nx̂, estim.nx̂, estim.nym, estim.He
     nX̂con = nx̂*(He+1)
-    notSolvedYet = (termination_status(optim) == OPTIMIZE_NOT_CALLED)
+    notSolvedYet = (JuMP.termination_status(optim) == JuMP.OPTIMIZE_NOT_CALLED)
     C = estim.C
     isnothing(X̂min)   && !isnothing(x̂min)   && (X̂min = repeat(x̂min, He+1))
     isnothing(X̂max)   && !isnothing(x̂max)   && (X̂max = repeat(x̂max, He+1))
@@ -507,8 +507,8 @@ function setconstraint!(
         A = con.A[con.i_b, :]
         b = con.b[con.i_b]
         Z̃var = optim[:Z̃var]
-        delete(optim, optim[:linconstraint])
-        unregister(optim, :linconstraint)
+        JuMP.delete(optim, optim[:linconstraint])
+        JuMP.unregister(optim, :linconstraint)
         @constraint(optim, linconstraint, A*Z̃var .≤ b)
         setnonlincon!(estim, model)
     else
@@ -578,22 +578,22 @@ setnonlincon!(::MovingHorizonEstimator, ::SimModel) = nothing
 function setnonlincon!(estim::MovingHorizonEstimator, ::NonLinModel)
     optim, con = estim.optim, estim.con
     Z̃var = optim[:Z̃var]
-    map(con -> delete(optim, con), all_nonlinear_constraints(optim))
+    map(con -> JuMP.delete(optim, con), JuMP.all_nonlinear_constraints(optim))
     for i in findall(.!isinf.(con.X̂min))
         f_sym = Symbol("g_X̂min_$(i)")
-        add_nonlinear_constraint(optim, :($(f_sym)($(Z̃var...)) <= 0))
+        JuMP.add_nonlinear_constraint(optim, :($(f_sym)($(Z̃var...)) <= 0))
     end
     for i in findall(.!isinf.(con.X̂max))
         f_sym = Symbol("g_X̂max_$(i)")
-        add_nonlinear_constraint(optim, :($(f_sym)($(Z̃var...)) <= 0))
+        JuMP.add_nonlinear_constraint(optim, :($(f_sym)($(Z̃var...)) <= 0))
     end
     for i in findall(.!isinf.(con.V̂min))
         f_sym = Symbol("g_V̂min_$(i)")
-        add_nonlinear_constraint(optim, :($(f_sym)($(Z̃var...)) <= 0))
+        JuMP.add_nonlinear_constraint(optim, :($(f_sym)($(Z̃var...)) <= 0))
     end
     for i in findall(.!isinf.(con.V̂max))
         f_sym = Symbol("g_V̂max_$(i)")
-        add_nonlinear_constraint(optim, :($(f_sym)($(Z̃var...)) <= 0))
+        JuMP.add_nonlinear_constraint(optim, :($(f_sym)($(Z̃var...)) <= 0))
     end
     return nothing
 end
@@ -967,7 +967,7 @@ function init_optimization!(
     estim::MovingHorizonEstimator, ::LinModel, optim::JuMP.GenericModel
 )
     nZ̃ = length(estim.Z̃)
-    set_silent(optim)
+    JuMP.set_silent(optim)
     limit_solve_time(estim.optim, estim.model.Ts)
     @variable(optim, Z̃var[1:nZ̃])
     A = estim.con.A[estim.con.i_b, :]
@@ -988,19 +988,19 @@ function init_optimization!(
     C, con = estim.C, estim.con
     nZ̃ = length(estim.Z̃)
     # --- variables and linear constraints ---
-    set_silent(optim)
+    JuMP.set_silent(optim)
     limit_solve_time(estim.optim, estim.model.Ts)
     @variable(optim, Z̃var[1:nZ̃])
     A = estim.con.A[con.i_b, :]
     b = estim.con.b[con.i_b]
     @constraint(optim, linconstraint, A*Z̃var .≤ b)
     # --- nonlinear optimization init ---
-    if !isinf(C) && solver_name(optim) == "Ipopt"
+    if !isinf(C) && JuMP.solver_name(optim) == "Ipopt"
         try
-            get_attribute(optim, "nlp_scaling_max_gradient")
+            JuMP.get_attribute(optim, "nlp_scaling_max_gradient")
         catch
             # default "nlp_scaling_max_gradient" to `10.0/C` if not already set:
-            set_attribute(optim, "nlp_scaling_max_gradient", 10.0/C)
+            JuMP.set_attribute(optim, "nlp_scaling_max_gradient", 10.0/C)
         end
     end
     He = estim.He