reduce allocation MovingHorizonEstimator

Author: franckgaga

.github/workflows/TagBot.yml

@@ -4,6 +4,22 @@ on:
     types:
       - created
   workflow_dispatch:
+    inputs:
+      lookback:
+        default: 3
+permissions:
+  actions: read
+  checks: read
+  contents: write
+  deployments: read
+  issues: read
+  discussions: read
+  packages: read
+  pages: read
+  pull-requests: read
+  repository-projects: read
+  security-events: read
+  statuses: read
 jobs:
   TagBot:
     if: github.event_name == 'workflow_dispatch' || github.actor == 'JuliaTagBot'
@@ -12,4 +28,6 @@ jobs:
       - uses: JuliaRegistries/TagBot@v1
         with:
           token: ${{ secrets.GITHUB_TOKEN }}
-          ssh: ${{ secrets.DOCUMENTER_KEY }}
+          # Edit the following line to reflect the actual name of the GitHub Secret containing your private key
+          ssh: ${{ secrets.DOCUMENTER_KEY }}
+          # ssh: ${{ secrets.NAME_OF_MY_SSH_PRIVATE_KEY_SECRET }}

README.md

@@ -105,6 +105,7 @@ for more detailed examples.
   - [x] quickly compare multiple optimizers
   - [x] nonlinear solvers relying on automatic differentiation (exact derivative)
 - [x] additional information about the optimum to ease troubleshooting
+- [x] implementation that carefully limits allocations for real-time applications
 
 ### State Estimation Features
 

src/controller/construct.jl

@@ -149,79 +149,79 @@ function setconstraint!(
     end
     if !isnothing(Umin)
         size(Umin)   == (nu*Hp,) || throw(ArgumentError("Umin size must be $((nu*Hp,))"))
-        con.Umin[:] = Umin
+        con.Umin .= Umin
     end
     if !isnothing(Umax)
         size(Umax)   == (nu*Hp,) || throw(ArgumentError("Umax size must be $((nu*Hp,))"))
-        con.Umax[:] = Umax
+        con.Umax .= Umax
     end
     if !isnothing(ΔUmin)
         size(ΔUmin)  == (nu*Hc,) || throw(ArgumentError("ΔUmin size must be $((nu*Hc,))"))
-        con.ΔŨmin[1:nu*Hc] = ΔUmin
+        con.ΔŨmin[1:nu*Hc] .= ΔUmin
     end
     if !isnothing(ΔUmax)
         size(ΔUmax)  == (nu*Hc,) || throw(ArgumentError("ΔUmax size must be $((nu*Hc,))"))
-        con.ΔŨmax[1:nu*Hc] = ΔUmax
+        con.ΔŨmax[1:nu*Hc] .= ΔUmax
     end
     if !isnothing(Ymin)
         size(Ymin)   == (ny*Hp,) || throw(ArgumentError("Ymin size must be $((ny*Hp,))"))
-        con.Ymin[:] = Ymin
+        con.Ymin .= Ymin
     end
     if !isnothing(Ymax)
         size(Ymax)   == (ny*Hp,) || throw(ArgumentError("Ymax size must be $((ny*Hp,))"))
-        con.Ymax[:] = Ymax
+        con.Ymax .= Ymax
     end
     if !isnothing(x̂min)
         size(x̂min)   == (nx̂,) || throw(ArgumentError("x̂min size must be $((nx̂,))"))
-        con.x̂min[:] = x̂min
+        con.x̂min .= x̂min
     end
     if !isnothing(x̂max)
         size(x̂max)   == (nx̂,) || throw(ArgumentError("x̂max size must be $((nx̂,))"))
-        con.x̂max[:] = x̂max
+        con.x̂max .= x̂max
     end
     if !isnothing(C_umin)
         size(C_umin) == (nu*Hp,) || throw(ArgumentError("C_umin size must be $((nu*Hp,))"))
         any(C_umin .< 0) && error("C_umin weights should be non-negative")
-        con.A_Umin[:, end] = -C_umin
+        con.A_Umin[:, end] .= -C_umin
     end
     if !isnothing(C_umax)
         size(C_umax) == (nu*Hp,) || throw(ArgumentError("C_umax size must be $((nu*Hp,))"))
         any(C_umax .< 0) && error("C_umax weights should be non-negative")
-        con.A_Umax[:, end] = -C_umax
+        con.A_Umax[:, end] .= -C_umax
     end
     if !isnothing(C_Δumin)
         size(C_Δumin) == (nu*Hc,) || throw(ArgumentError("C_Δumin size must be $((nu*Hc,))"))
         any(C_Δumin .< 0) && error("C_Δumin weights should be non-negative")
-        con.A_ΔŨmin[1:end-1, end] = -C_Δumin 
+        con.A_ΔŨmin[1:end-1, end] .= -C_Δumin 
     end
     if !isnothing(C_Δumax)
         size(C_Δumax) == (nu*Hc,) || throw(ArgumentError("C_Δumax size must be $((nu*Hc,))"))
         any(C_Δumax .< 0) && error("C_Δumax weights should be non-negative")
-        con.A_ΔŨmax[1:end-1, end] = -C_Δumax
+        con.A_ΔŨmax[1:end-1, end] .= -C_Δumax
     end
     if !isnothing(C_ymin)
         size(C_ymin) == (ny*Hp,) || throw(ArgumentError("C_ymin size must be $((ny*Hp,))"))
         any(C_ymin .< 0) && error("C_ymin weights should be non-negative")
-        con.C_ymin[:] = C_ymin
-        size(con.A_Ymin, 1) ≠ 0 && (con.A_Ymin[:, end] = -con.C_ymin) # for LinModel
+        con.C_ymin .= C_ymin
+        size(con.A_Ymin, 1) ≠ 0 && (con.A_Ymin[:, end] .= -con.C_ymin) # for LinModel
     end
     if !isnothing(C_ymax)
         size(C_ymax) == (ny*Hp,) || throw(ArgumentError("C_ymax size must be $((ny*Hp,))"))
         any(C_ymax .< 0) && error("C_ymax weights should be non-negative")
-        con.C_ymax[:] = C_ymax
-        size(con.A_Ymax, 1) ≠ 0 && (con.A_Ymax[:, end] = -con.C_ymax) # for LinModel
+        con.C_ymax .= C_ymax
+        size(con.A_Ymax, 1) ≠ 0 && (con.A_Ymax[:, end] .= -con.C_ymax) # for LinModel
     end
     if !isnothing(c_x̂min)
         size(c_x̂min) == (nx̂,) || throw(ArgumentError("c_x̂min size must be $((nx̂,))"))
         any(c_x̂min .< 0) && error("c_x̂min weights should be non-negative")
-        con.c_x̂min[:] = c_x̂min
-        size(con.A_x̂min, 1) ≠ 0 && (con.A_x̂min[:, end] = -con.c_x̂min) # for LinModel
+        con.c_x̂min .= c_x̂min
+        size(con.A_x̂min, 1) ≠ 0 && (con.A_x̂min[:, end] .= -con.c_x̂min) # for LinModel
     end
     if !isnothing(c_x̂max)
         size(c_x̂max) == (nx̂,) || throw(ArgumentError("c_x̂max size must be $((nx̂,))"))
         any(c_x̂max .< 0) && error("c_x̂max weights should be non-negative")
-        con.c_x̂max[:] = c_x̂max
-        size(con.A_x̂max, 1) ≠ 0 && (con.A_x̂max[:, end] = -con.c_x̂max) # for LinModel
+        con.c_x̂max .= c_x̂max
+        size(con.A_x̂max, 1) ≠ 0 && (con.A_x̂max[:, end] .= -con.c_x̂max) # for LinModel
     end
     i_Umin,  i_Umax  = .!isinf.(con.Umin),  .!isinf.(con.Umax)
     i_ΔŨmin, i_ΔŨmax = .!isinf.(con.ΔŨmin), .!isinf.(con.ΔŨmin)

src/controller/execute.jl

@@ -225,15 +225,15 @@ function predictstoch!(
 ) where {NT<:Real}
     isnothing(ym) && error("Predictive controllers with InternalModel need the measured "*
                            "outputs ym in keyword argument to compute control actions u")
-    Ŷop = mpc.Ŷop
+    Ŷop, ny, yop = mpc.Ŷop, estim.model.ny, estim.model.yop
     ŷd = h(estim.model, estim.x̂d, d - estim.model.dop) .+ estim.model.yop 
     ŷs = zeros(NT, estim.model.ny)
     ŷs[estim.i_ym] .= @views ym .- ŷd[estim.i_ym]  # ŷs=0 for unmeasured outputs
     Ŷop_LHS = similar(Ŷop)
     Ŷop  .= mul!(Ŷop_LHS, mpc.Ks, estim.x̂s)
     Ŷop .+= mul!(Ŷop_LHS, mpc.Ps, ŷs)
     for j=1:mpc.Hp
-        Ŷop[(1 + estim.model.ny*(j-1)):(estim.model.ny*j)] .+= estim.model.yop
+        Ŷop[(1 + ny*(j-1)):(ny*j)] .+= yop
     end
     return nothing
 end
@@ -249,6 +249,8 @@ Also init ``\mathbf{f_x̂} = \mathbf{g_x̂ d}(k) + \mathbf{j_x̂ D̂} + \mathbf{
 vector for the terminal constraints, see [`init_predmat`](@ref).
 """
 function linconstraint!(mpc::PredictiveController, model::LinModel)
+    nU, nΔŨ, nY = length(mpc.con.Umin), length(mpc.con.ΔŨmin), length(mpc.con.Ymin)
+    nx̂ = mpc.estim.nx̂
     fx̂ = mpc.con.fx̂
     fx̂_LHS = similar(fx̂)
     fx̂  .= mul!(fx̂_LHS, mpc.con.kx̂, mpc.estim.x̂)
@@ -257,28 +259,37 @@ function linconstraint!(mpc::PredictiveController, model::LinModel)
         fx̂ .+= mul!(fx̂_LHS, mpc.con.gx̂, mpc.d0)
         fx̂ .+= mul!(fx̂_LHS, mpc.con.jx̂, mpc.D̂0)
     end
-    mpc.con.b .= [
-        -mpc.con.Umin + mpc.T_lastu
-        +mpc.con.Umax - mpc.T_lastu
-        -mpc.con.ΔŨmin
-        +mpc.con.ΔŨmax 
-        -mpc.con.Ymin + mpc.F
-        +mpc.con.Ymax - mpc.F
-        -mpc.con.x̂min + mpc.con.fx̂
-        +mpc.con.x̂max - mpc.con.fx̂
-    ]
+    n = 0
+    mpc.con.b[(n+1):(n+nU)]  .= @. -mpc.con.Umin + mpc.T_lastu
+    n += nU
+    mpc.con.b[(n+1):(n+nU)]  .= @. +mpc.con.Umax - mpc.T_lastu
+    n += nU
+    mpc.con.b[(n+1):(n+nΔŨ)] .= @. -mpc.con.ΔŨmin
+    n += nΔŨ
+    mpc.con.b[(n+1):(n+nΔŨ)] .= @. +mpc.con.ΔŨmax
+    n += nΔŨ
+    mpc.con.b[(n+1):(n+nY)]  .= @. -mpc.con.Ymin + mpc.F
+    n += nY
+    mpc.con.b[(n+1):(n+nY)]  .= @. +mpc.con.Ymax - mpc.F
+    n += nY
+    mpc.con.b[(n+1):(n+nx̂)]  .= @. -mpc.con.x̂min + fx̂
+    n += nx̂
+    mpc.con.b[(n+1):(n+nx̂)]  .= @. +mpc.con.x̂max - fx̂
     lincon = mpc.optim[:linconstraint]
     set_normalized_rhs.(lincon, mpc.con.b[mpc.con.i_b])
 end
 
 "Set `b` excluding predicted output constraints when `model` is not a [`LinModel`](@ref)."
-function linconstraint!(mpc::PredictiveController, model::SimModel)
-    mpc.con.b .= [
-        -mpc.con.Umin + mpc.T_lastu
-        +mpc.con.Umax - mpc.T_lastu
-        -mpc.con.ΔŨmin
-        +mpc.con.ΔŨmax 
-    ]
+function linconstraint!(mpc::PredictiveController, ::SimModel)
+    nU, nΔŨ = length(mpc.con.Umin), length(mpc.con.ΔŨmin)
+    n = 0
+    mpc.con.b[(n+1):(n+nU)]  .= @. -mpc.con.Umin + mpc.T_lastu
+    n += nU
+    mpc.con.b[(n+1):(n+nU)]  .= @. +mpc.con.Umax - mpc.T_lastu
+    n += nU
+    mpc.con.b[(n+1):(n+nΔŨ)] .= @. -mpc.con.ΔŨmin
+    n += nΔŨ
+    mpc.con.b[(n+1):(n+nΔŨ)] .= @. +mpc.con.ΔŨmax
     lincon = mpc.optim[:linconstraint]
     set_normalized_rhs.(lincon, mpc.con.b[mpc.con.i_b])
 end
@@ -442,7 +453,7 @@ function optim_objective!(mpc::PredictiveController{NT}) where {NT<:Real}
         end
         @debug solution_summary(optim, verbose=true)
     end
-    mpc.ΔŨ[:] = iserror(optim) ? ΔŨlast : ΔŨcurr
+    mpc.ΔŨ .= iserror(optim) ? ΔŨlast : ΔŨcurr
     return mpc.ΔŨ
 end
 

src/estimator/kalman.jl

@@ -752,7 +752,7 @@ automatically computes the Jacobians:
 ```math
 \begin{aligned}
     \mathbf{F̂}(k) &= \left. \frac{∂\mathbf{f̂}(\mathbf{x̂}, \mathbf{u}, \mathbf{d})}{∂\mathbf{x̂}} \right|_{\mathbf{x̂ = x̂}_{k-1}(k),\, \mathbf{u = u}(k),\, \mathbf{d = d}(k)}  \\
-    \mathbf{Ĥ}(k) &= \left. \frac{∂\mathbf{ĥ}(\mathbf{x̂}, \mathbf{d})}{∂\mathbf{x̂}}             \right|_{\mathbf{x = x̂}_{k-1}(k),\, \mathbf{d = d}(k)}
+    \mathbf{Ĥ}(k) &= \left. \frac{∂\mathbf{ĥ}(\mathbf{x̂}, \mathbf{d})}{∂\mathbf{x̂}}             \right|_{\mathbf{x̂ = x̂}_{k-1}(k),\, \mathbf{d = d}(k)}
 \end{aligned}
 ```
 The matrix ``\mathbf{Ĥ^m}`` is the rows of ``\mathbf{Ĥ}`` that are measured outputs.

src/estimator/mhe/construct.jl

@@ -418,63 +418,65 @@ function setconstraint!(
     if !isnothing(X̂min)
         size(X̂min) == (nX̂con,) || throw(ArgumentError("X̂min size must be $((nX̂con,))"))
         con.x̃min[end-nx̂+1:end] = X̂min[1:nx̂] # if C is finite : x̃ = [ϵ; x̂]
-        con.X̂min[:] = X̂min[nx̂+1:end]
+        con.X̂min .= X̂min[nx̂+1:end]
     end
     if !isnothing(X̂max)
         size(X̂max) == (nX̂con,) || throw(ArgumentError("X̂max size must be $((nX̂con,))"))
         con.x̃max[end-nx̂+1:end] = X̂max[1:nx̂] # if C is finite : x̃ = [ϵ; x̂]
-        con.X̂max[:] = X̂max[nx̂+1:end]
+        con.X̂max .= X̂max[nx̂+1:end]
     end
     if !isnothing(Ŵmin)
         size(Ŵmin) == (nŵ*He,) || throw(ArgumentError("Ŵmin size must be $((nŵ*He,))"))
-        con.Ŵmin[:] = Ŵmin
+        con.Ŵmin .= Ŵmin
     end
     if !isnothing(Ŵmax)
         size(Ŵmax) == (nŵ*He,) || throw(ArgumentError("Ŵmax size must be $((nŵ*He,))"))
-        con.Ŵmax[:] = Ŵmax
+        con.Ŵmax .= Ŵmax
     end
     if !isnothing(V̂min)
         size(V̂min) == (nym*He,) || throw(ArgumentError("V̂min size must be $((nym*He,))"))
-        con.V̂min[:] = V̂min
+        con.V̂min .= V̂min
     end
     if !isnothing(V̂max)
         size(V̂max) == (nym*He,) || throw(ArgumentError("V̂max size must be $((nym*He,))"))
-        con.V̂max[:] = V̂max
+        con.V̂max .= V̂max
     end
     if !isnothing(C_x̂min)
         size(C_x̂min) == (nX̂con,) || throw(ArgumentError("C_x̂min size must be $((nX̂con,))"))
         any(C_x̂min .< 0) && error("C_x̂min weights should be non-negative")
-        con.A_x̃min[end-nx̂+1:end, end] = -C_x̂min[1:nx̂] # if C is finite : x̃ = [ϵ; x̂]
-        con.C_x̂min[:] = C_x̂min[nx̂+1:end]
+        # if C is finite : x̃ = [ϵ; x̂] 
+        con.A_x̃min[end-nx̂+1:end, end] .= @views -C_x̂min[1:nx̂] 
+        con.C_x̂min .= @views C_x̂min[nx̂+1:end]
         size(con.A_X̂min, 1) ≠ 0 && (con.A_X̂min[:, end] = -con.C_x̂min) # for LinModel
     end
     if !isnothing(C_x̂max)
         size(C_x̂max) == (nX̂con,) || throw(ArgumentError("C_x̂max size must be $((nX̂con,))"))
         any(C_x̂max .< 0) && error("C_x̂max weights should be non-negative")
-        con.A_x̃max[end-nx̂+1:end, end] = -C_x̂max[1:nx̂] # if C is finite : x̃ = [ϵ; x̂]
-        con.C_x̂max[:] = C_x̂max[nx̂+1:end]
+        # if C is finite : x̃ = [ϵ; x̂] :
+        con.A_x̃max[end-nx̂+1:end, end] .= @views -C_x̂max[1:nx̂]
+        con.C_x̂max .= @views C_x̂max[nx̂+1:end]
         size(con.A_X̂max, 1) ≠ 0 && (con.A_X̂max[:, end] = -con.C_x̂max) # for LinModel
     end
     if !isnothing(C_ŵmin)
         size(C_ŵmin) == (nŵ*He,) || throw(ArgumentError("C_ŵmin size must be $((nŵ*He,))"))
         any(C_ŵmin .< 0) && error("C_ŵmin weights should be non-negative")
-        con.A_Ŵmin[:, end] = -C_ŵmin
+        con.A_Ŵmin[:, end] .= -C_ŵmin
     end
     if !isnothing(C_ŵmax)
         size(C_ŵmax) == (nŵ*He,) || throw(ArgumentError("C_ŵmax size must be $((nŵ*He,))"))
         any(C_ŵmax .< 0) && error("C_ŵmax weights should be non-negative")
-        con.A_Ŵmax[:, end] = -C_ŵmax
+        con.A_Ŵmax[:, end] .= -C_ŵmax
     end
     if !isnothing(C_v̂min)
         size(C_v̂min) == (nym*He,) || throw(ArgumentError("C_v̂min size must be $((nym*He,))"))
         any(C_v̂min .< 0) && error("C_v̂min weights should be non-negative")
-        con.C_v̂min[:] = C_v̂min
+        con.C_v̂min .= C_v̂min
         size(con.A_V̂min, 1) ≠ 0 && (con.A_V̂min[:, end] = -con.C_v̂min) # for LinModel
     end
     if !isnothing(C_v̂max)
         size(C_v̂max) == (nym*He,) || throw(ArgumentError("C_v̂max size must be $((nym*He,))"))
         any(C_v̂max .< 0) && error("C_v̂max weights should be non-negative")
-        con.C_v̂max[:] = C_v̂max
+        con.C_v̂max .= C_v̂max
         size(con.A_V̂max, 1) ≠ 0 && (con.A_V̂max[:, end] = -con.C_v̂max) # for LinModel
     end
     i_x̃min, i_x̃max  = .!isinf.(con.x̃min)  , .!isinf.(con.x̃max)
@@ -989,40 +991,46 @@ function init_optimization!(
     He = estim.He
     nV̂, nX̂, ng = He*estim.nym, He*estim.nx̂, length(con.i_g)
     # see init_optimization!(mpc::NonLinMPC, optim) for details on the inspiration
-    Jfunc, gfunc = let estim=estim, model=model, nZ̃=nZ̃ , nV̂=nV̂, nX̂=nX̂, ng=ng
+    Jfunc, gfunc = let estim=estim, model=model, nZ̃=nZ̃ , nV̂=nV̂, nX̂=nX̂, ng=ng, nx̂=estim.nx̂
         last_Z̃tup_float, last_Z̃tup_dual = nothing, nothing
         V̂_cache::DiffCache{Vector{JNT}, Vector{JNT}} = DiffCache(zeros(JNT, nV̂), nZ̃ + 3)
         g_cache::DiffCache{Vector{JNT}, Vector{JNT}} = DiffCache(zeros(JNT, ng), nZ̃ + 3)
         X̂_cache::DiffCache{Vector{JNT}, Vector{JNT}} = DiffCache(zeros(JNT, nX̂), nZ̃ + 3)
+        x̄_cache::DiffCache{Vector{JNT}, Vector{JNT}} = DiffCache(zeros(JNT, nx̂), nZ̃ + 3)
         function Jfunc(Z̃tup::JNT...)
-            V̂ = get_tmp(V̂_cache, Z̃tup[1])
+            Z̃1 = Z̃tup[begin]
+            V̂ = get_tmp(V̂_cache, Z̃1)
             Z̃ = collect(Z̃tup)
             if Z̃tup !== last_Z̃tup_float
-                g = get_tmp(g_cache, Z̃tup[1])
-                X̂ = get_tmp(X̂_cache, Z̃tup[1])
+                g = get_tmp(g_cache, Z̃1)
+                X̂ = get_tmp(X̂_cache, Z̃1)
                 V̂, X̂ = predict!(V̂, X̂, estim, model, Z̃)
                 g = con_nonlinprog!(g, estim, model, X̂, V̂, Z̃)
                 last_Z̃tup_float = Z̃tup
             end
-            return obj_nonlinprog(estim, model, V̂, Z̃)
+            x̄ = get_tmp(x̄_cache, Z̃1)
+            return obj_nonlinprog!(x̄, estim, model, V̂, Z̃)
         end
         function Jfunc(Z̃tup::ForwardDiff.Dual...)
-            V̂ = get_tmp(V̂_cache, Z̃tup[1])
+            Z̃1 = Z̃tup[begin]
+            V̂ = get_tmp(V̂_cache, Z̃1)
             Z̃ = collect(Z̃tup)
             if Z̃tup !== last_Z̃tup_dual
-                g = get_tmp(g_cache, Z̃tup[1])
-                X̂ = get_tmp(X̂_cache, Z̃tup[1])
+                g = get_tmp(g_cache, Z̃1)
+                X̂ = get_tmp(X̂_cache, Z̃1)
                 V̂, X̂ = predict!(V̂, X̂, estim, model, Z̃)
                 g = con_nonlinprog!(g, estim, model, X̂, V̂, Z̃)
                 last_Z̃tup_dual = Z̃tup
             end
-            return obj_nonlinprog(estim, model, V̂, Z̃)
+            x̄ = get_tmp(x̄_cache, Z̃1)
+            return obj_nonlinprog!(x̄, estim, model, V̂, Z̃)
         end
         function gfunc_i(i, Z̃tup::NTuple{N, JNT}) where N
-            g = get_tmp(g_cache, Z̃tup[1])
+            Z̃1 = Z̃tup[begin]
+            g = get_tmp(g_cache, Z̃1)
             if Z̃tup !== last_Z̃tup_float
-                V̂ = get_tmp(V̂_cache, Z̃tup[1])
-                X̂ = get_tmp(X̂_cache, Z̃tup[1])
+                V̂ = get_tmp(V̂_cache, Z̃1)
+                X̂ = get_tmp(X̂_cache, Z̃1)
                 Z̃ = collect(Z̃tup)
                 V̂, X̂ = predict!(V̂, X̂, estim, model, Z̃)
                 g = con_nonlinprog!(g, estim, model, X̂, V̂, Z̃)
@@ -1031,10 +1039,11 @@ function init_optimization!(
             return g[i]
         end 
         function gfunc_i(i, Z̃tup::NTuple{N, ForwardDiff.Dual}) where N
-            g = get_tmp(g_cache, Z̃tup[1])
+            Z̃1 = Z̃tup[begin]
+            g = get_tmp(g_cache, Z̃1)
             if Z̃tup !== last_Z̃tup_dual
-                V̂ = get_tmp(V̂_cache, Z̃tup[1])
-                X̂ = get_tmp(X̂_cache, Z̃tup[1])
+                V̂ = get_tmp(V̂_cache, Z̃1)
+                X̂ = get_tmp(X̂_cache, Z̃1)
                 Z̃ = collect(Z̃tup)
                 V̂, X̂ = predict!(V̂, X̂, estim, model, Z̃)
                 g = con_nonlinprog!(g, estim, model, X̂, V̂, Z̃)