# Help solving: Adding a few equations to ODE slows JuMP >100x

I’m using JuMP with Ipopt and AppleAccelerate.jl

Following the rocket control tutorial, I coded a simplified set of differential equations where I optimize a control function. My simplified code takes <2s to run:

	# dynamical variables
u = @JuMP.variables(model,begin
0 <= φ_R[1:n] <= 0.45
0 <= c_pc[1:n]
0 <= M[1:n]
end);

# optimization knob
@variable(model, 0 <= α_R[1:n] <= 0.45 )

# helper functions
@expression(model, γ[j = 1:n], 10 * c_pc[j]/(c_pc[j] + 0.03))
@expression(model, λ[j = 1:n], γ[j] * φ_R[j])

du = Vector{Vector{NonlinearExpr}}(undef, length(u))
# growth dynamics
du[1] = @expression(model, δφ_R[j = 1:n], (α_R[j] - φ_R[j]) * λ[j])
du[2] = @expression(model, δc_pc[j =1:n], νmax*(0.45 - φ_R[j]) - (1 + c_pc[j]) * λ[j] )
du[3] = @expression(model, δM[j = 1:n], M[j] * λ[j])


(and with trapezoidal integration)

the expanded model takes ~250s to run:

	# dynamical variables
u = @JuMP.variables(model,begin
0 <= φ_R[1:n] <= 0.45
0 <= c_pc[1:n]
0 <= C_bnd[1:n]
0 <= C_in[1:n]
0 <= M[1:n]
end);

# optimization knob
@variable(model, 0 <= α_R[1:n] <= 0.45 )

# helper functions
C_max = 1e-5
@expression(model, C_out[j = 1:n], C_max)
@expression(model, C_Rf[j = 1:n], φ_R[j] * 6.4e-5 - C_bnd[j])
@expression(model, φ_Rf[j = 1:n], φ_R[j] - C_bnd[j] / 6.4e-5)
@expression(model, φ_M[j = 1:n], 0.45 - φ_R[j])
@expression(model, γ[j = 1:n], 10 * c_pc[j]/(c_pc[j] + 0.03))
@expression(model, λ[j = 1:n], γ[j] * φ_Rf[j])

du = Vector{Vector{NonlinearExpr}}(undef, length(u))
# growth dynamics
k_on =  1.1 # 4e5
du[1] = @expression(model, δφ_R[j = 1:n], (α_R[j] - φ_R[j]) * λ[j])
du[2] = @expression(model, δc_pc[j =1:n], 5*(0.45 - φ_R[j]) - (1 + c_pc[j]) * λ[j] )
du[3] = @expression(model, δC_bnd[j = 1:n], k_on*C_in[j]*C_Rf[j] - 0.4*C_bnd[j] - λ[j]*C_bnd[j])
du[4] = @expression(model, δC_in[j = 1:n],
1.5*(C_out[j] - C_in[j]) + 0.4*C_bnd[j] - k_on*C_Rf[j]*C_in[j] - λ[j]*C_in[j]
)
du[5] = @expression(model, δM[j = 1:n], M[j] * λ[j])


There is an even more expanded model with another two equations that I haven’t been patient enough to get to converge at all. I previously solved these equations in the DifferentialEquations package, so barring typos, I’m reasonably confident that the equations are not the problem.

Any suggestions on how to speed up finding the optimal control law would be much appreciated!

Can you provide a reproducible example?

The actual numerics functions

begin
abstract type AbstractIntegrationRule end
struct TrapezoidalIntegrationRule <: AbstractIntegrationRule end

struct JumpConfig{T<:AbstractIntegrationRule}
n::Int   # mesh points
Δt::Float64  # time step
integration_rule::T
end

function integration_JuMP!(model, u, du, jc, ::TrapezoidalIntegrationRule)
@variable(model, δt[1:jc.n] == jc.Δt)

for t in 2:jc.n
s = t - 1
for k in 1:length(u)
@constraint(model, u[k][t] == u[k][s] + 0.5 * δt[t] * (du[k][s] + du[k][t]))
end
end
end

function fast_sim(
u0, jc::JumpConfig
)
n = jc.n
model = Model(Ipopt.Optimizer)

# dynamical variables
u = @JuMP.variables(model,begin
0 <= φ_R[1:n] <= 0.45
0 <= c_pc[1:n]
0 <= M[1:n]
end);

# optimization knob
@variable(model, 0 <= α_R[1:n] <= 0.45 )

# helper functions
@expression(model, γ[j = 1:n], 10 * c_pc[j]/(c_pc[j] + 0.03))
@expression(model, λ[j = 1:n], γ[j] * φ_R[j])

du = Vector{Vector{NonlinearExpr}}(undef, length(u))
# growth dynamics
du[1] = @expression(model, δφ_R[j = 1:n], (α_R[j] - φ_R[j]) * λ[j])
du[2] = @expression(model, δc_pc[j =1:n], 5*(0.45 - φ_R[j]) - (1 + c_pc[j]) * λ[j] )
du[3] = @expression(model, δM[j = 1:n], M[j] * λ[j])

# integration
integration_JuMP!(model, u, du, jc, TrapezoidalIntegrationRule())

@objective(model, Max, M[n])
set_silent(model)

# one set of initial conditions
for i in 1:length(u)
fix(u[i][1], u0[i]; force = true)
end

optimize!(model)
Ls = value.(λ)
Rs = value.(φ_R)
aRs = value.(α_R)
cpcs = value.(c_pc)

return Rs, cpcs, aRs, Ls
end

function slow_sim(
u0, jc::JumpConfig; Kd = 1e-6, k_on = 4e6, C_max = 1e-6
)
n = jc.n
model = Model(Ipopt.Optimizer)

# dynamical variables
u = @JuMP.variables(model,begin
0 <= φ_R[1:n] <= 0.45
0 <= c_pc[1:n]
0 <= C_bnd[1:n]
0 <= C_in[1:n]
0 <= M[1:n]
end);

# optimization knob
@variable(model, 0 <= α_R[1:n] <= 0.45 )

# helper functions
@expression(model, C_out[j = 1:n], C_max)
@expression(model, C_Rf[j = 1:n], φ_R[j] * 6.4e-5 - C_bnd[j])

@expression(model, φ_Rf[j = 1:n], φ_R[j] - C_bnd[j] / 6.4e-5)
@constraint(model, φ_Rf .>= 0)

@expression(model, φ_M[j = 1:n], 0.45 - φ_R[j])
@expression(model, γ[j = 1:n], 10 * c_pc[j]/(c_pc[j] + 0.03))
@expression(model, λ[j = 1:n], γ[j] * φ_Rf[j])

du = Vector{Vector{NonlinearExpr}}(undef, length(u))
# growth dynamics
k_off = k_on * Kd # 0.4  # 40
du[1] = @expression(model, δφ_R[j = 1:n], (α_R[j] - φ_R[j]) * λ[j])
# 1/s in mass fractions -> multiply by ρ to get conc.
du[2] = @expression(model, δc_pc[j =1:n], 5*(0.45 - φ_R[j]) - (1 + c_pc[j]) * λ[j] )
du[3] = @expression(model, δC_bnd[j = 1:n], k_on*C_in[j]*C_Rf[j] - k_off*C_bnd[j] - λ[j]*C_bnd[j])
du[4] = @expression(model, δC_in[j = 1:n],
1.5*(C_out[j] - C_in[j]) + k_off*C_bnd[j] - k_on*C_Rf[j]*C_in[j] - λ[j]*C_in[j]
)
du[5] = @expression(model, δM[j = 1:n], M[j] * λ[j])

# integration
integration_JuMP!(model, u, du, jc, TrapezoidalIntegrationRule())

@objective(model, Max, M[n])
set_silent(model)

for i in 1:length(u)
fix(u[i][1], u0[i]; force = true)
end

optimize!(model)
Ls = value.(λ)
Rs = value.(φ_R)
a_Rs = value.(α_R)
cpcs = value.(c_pc)
C_bnds = value.(C_bnd)
C_ins = value.(C_in)

return Rs, a_Rs, cpcs, C_bnds, C_ins, Ls
end
end


Example code to run it:

begin
u0_fast = [0.1, 0.05, 0.001]

steps = 1000
δt = 1/250
jc = JumpConfig(steps, δt, TrapezoidalIntegrationRule())
fast_soln = fast_sim(u0_fast, jc)
end


It takes ~2s to run the above in a Pluto NB ^

u0_slow = [0.1, 0.05, 0.0, 0.0, 0.001]
slow_soln = slow_sim(u0_slow, jc)


and ~30s to run ^ the slow sim.

I was able to speed it up by adding a constraint, \lambda \geq 0, but still looking for it to be faster.

The speed is significantly affected by the magnitudes of k_on and C_max; I have them set lower than I want them to get the simulation to run. Something about the significant difference in the parameter magnitudes is significantly slowing down the optimization.

Here are also some plotting utilities for plotting the slow solution, in case it is useful. The results from the above optimizations can be plotted like plot_soln(slow_soln, jc).

My goal is to figure out how to get the slow_sim faster, fast_sim was just an intermediate step that I am sharing for background.

const green = colorant"#04C13C"
const purp = colorant"#BA79BE"
const pink = colorant"#FA007B"
const blue = colorant"#05A8C9"
const brown = colorant"#8B4513"
const orange = colorant"#FFA500"

const thick_line = 3.5
const thin_line = 2
const extra_thin_line = 1
const opt_color = :gray

function plot_soln(soln, jc::JumpConfig)
Rs = soln[1]
aRs = soln[2]
cpcs = soln[3]
normed_cpcs = cpcs ./ 0.03
C_bnds = soln[4]
C_ins = soln[5]
Ls = soln[6]

ts = cumsum([0; repeat([jc.Δt], jc.n)])[1:end-1]

plt_growth =  λ_plot(ts, Ls)
plt_pfractions = φR_plot(ts, Rs, aRs)
plt_pcursors = precursor_plot(ts, normed_cpcs, xlab=true)
plt_ABX = ABX_plot(ts, C_ins, C_bnds, xlab=true)

p = plot(
plt_pfractions, plt_growth, plt_pcursors, plt_ABX;
layout = grid(2,2), size = (700,700), legend=:right
)
return p
end

function precursor_plot(ts, normed_c_pcs; subset_idxs=nothing, xlab=false)
if xlab
xlab = "Time (hrs)"
else
xlab = ""
end
plt = plot_series(
ts, normed_c_pcs, blue;
title = "Precursors",
xlab = xlab,
ylab = "Concentration (Km)",
subset_idxs = subset_idxs,
)
return plt
end

function φR_plot(ts, φ_Rs, α_Rs; subset_idxs=nothing, xlab=false)
if xlab
xlab = "Time (hrs)"
else
xlab = ""
end
plt = plot_series(
ts, φ_Rs, purp;
title = "Protein Fractions",
xlab = xlab,
ylab = "Mass Fraction",
subset_idxs = subset_idxs,
)
plot_series!(
plt, ts, α_Rs, pink;
subset_idxs = subset_idxs, lw = extra_thin_line, ls = :dashdot,
)
return plt
end

function λ_plot(ts, λs; subset_idxs=nothing, xlab=false)
if xlab
xlab = "Time (hrs)"
else
xlab = ""
end
plt = plot_series(
ts, λs, green;
title = "Growth Rate",
xlab = xlab,
ylab = "Rate (1/hr)",
subset_idxs = subset_idxs,
)
return plt
end

function ABX_plot(ts, C_ins, C_bnds; subset_idxs=nothing, xlab=false)
if xlab
xlab = "Time (hrs)"
else
xlab = ""
end
plt = plot_series(
ts, C_bnds, orange;
title = "ABX",
xlab = xlab,
ylab = "Concentration ([M])",
subset_idxs = subset_idxs,
label = "Bound to Ribosome"
)
plot_series!(
plt, ts, C_ins, brown;
subset_idxs = subset_idxs,
label = "Intracellular",
)

return plt
end

function plot_series(
ts, xs, color::RGB;
title, xlab, ylab,
subset_idxs=nothing, lw = thick_line,
ls = :solid, label = nothing
)
xs, n_plots = subset_series(xs, subset_idxs)
colors = set_colors(n_plots, color)

plt = plot(title = title, xlabel = xlab, ylabel = ylab)
for i in 1:n_plots
if i == n_plots
lab = label
else
lab = nothing
end
plot!(
ts,
xs[:,i];
label = lab,
linecolor = colors[i],
linewidth = lw, linestyle = ls,
)
end
return plt
end

function plot_series!(
plt, ts, xs, color::RGB;
subset_idxs=nothing, lw = thick_line,
ls = :solid, label = nothing
)
xs, n_plots = subset_series(xs, subset_idxs)
colors = set_colors(n_plots, color)

for i in 1:n_plots
if i == n_plots
lab = label
else
lab = nothing
end
plot!(
plt,
ts,
xs[:,i];
label = lab,
linecolor = colors[i],
linewidth = lw, linestyle = ls,
)
end
return plt
end

function set_colors(n_plots, base)
if n_plots == 1
return [base]
else
return gen_color_range(n_plots, base)
end
end

function subset_series(series, subset_idxs)
if series isa Matrix
if subset_idxs != nothing
series = series[:, subset_idxs]
end
n_plots = size(series)[2]
elseif series isa Vector
n_plots = 1
else
error("series must be a matrix or vector")
end
return series, n_plots
end

function gen_color_range(n, end_color, start_color=colorant"cornsilk3")
if n == 1
return [end_color]
end
# rescale
xs = (0:n-1)/(5/4*(n-1)) .+ 0.2
cs = ColorScheme(range(start_color, end_color))
return get(cs, xs)
end


Can you provide the inputs for fast_sim and slow_sim?

Have you looked at the Ipopt logs? What part is taking the time?

how do I look at the Ipopt logs?

my intution is that it’s taking a long time due to the fact that the added DEs have terms that have very large and very small magnitudes, that also feedback onto the rest of the dynamics

how do I look at the Ipopt logs?

Comment out set_silent(model)

I get ERROR: LoadError: UndefVarError: νmax not defined.

I fixed it in the above code, but you can replace it with a 5. I originally was passing in a struct but removed it to minimize extra code

The Ipopt log gives me:

(I cut out some lines in the middle of the objective changing bc it exceeded the max character count)

This is Ipopt version 3.14.14, running with linear solver MUMPS 5.6.2.

Number of nonzeros in equality constraint Jacobian...:    37944
Number of nonzeros in inequality constraint Jacobian.:     1998
Number of nonzeros in Lagrangian Hessian.............:   164721

Total number of variables............................:     5995
variables with only lower bounds:     3996
variables with lower and upper bounds:     1999
variables with only upper bounds:        0
Total number of equality constraints.................:     4995
Total number of inequality constraints...............:     1000
inequality constraints with only lower bounds:     1000
inequality constraints with lower and upper bounds:        0
inequality constraints with only upper bounds:        0

iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
0  9.9999900e-03 1.56e+02 2.00e+00  -1.0 0.00e+00    -  0.00e+00 0.00e+00   0
1r 9.9999900e-03 1.56e+02 9.99e+02   0.2 0.00e+00    -  0.00e+00 3.33e-07R  4
2r 1.9263008e-02 6.12e+01 9.98e+02   0.2 1.65e+03    -  1.22e-03 9.70e-04f  1
3  1.9262650e-02 6.12e+01 3.20e+03  -1.0 6.13e+01    -  9.73e-05 1.69e-04f  1
4  1.9262614e-02 6.12e+01 3.06e+03  -1.0 3.36e+02    -  4.14e-05 1.79e-05h  1
5  1.9262503e-02 6.12e+01 2.87e+03  -1.0 5.90e+02    -  3.38e-05 5.45e-05h  1
6  1.9262502e-02 6.12e+01 2.87e+03  -1.0 6.18e+02    -  1.23e-04 2.37e-07h  6
7  1.9262502e-02 6.12e+01 2.87e+03  -1.0 5.82e+02    -  1.32e-04 2.87e-07h  6
8  1.9262465e-02 6.12e+01 2.87e+03  -1.0 5.45e+02    -  6.34e-05 1.85e-05h  1
9  1.9262456e-02 6.12e+01 2.87e+03  -1.0 4.80e+03    -  2.70e-06 4.24e-06f  2
iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
10  1.9262456e-02 6.12e+01 2.99e+03  -1.0 4.90e+03    -  6.86e-06 2.63e-07h  5
11  1.9262452e-02 6.12e+01 3.29e+03  -1.0 4.88e+03    -  4.94e-06 1.60e-06f  3
12  1.9262451e-02 6.12e+01 3.83e+03  -1.0 4.85e+03    -  6.00e-06 6.11e-07h  4
13  1.9262444e-02 6.12e+01 3.48e+03  -1.0 4.81e+03    -  5.26e-07 3.86e-06f  2
14  1.9262436e-02 6.12e+01 3.30e+03  -1.0 4.86e+03    -  2.14e-06 3.72e-06f  2
15r 1.9262436e-02 6.12e+01 1.00e+03  -0.4 0.00e+00    -  0.00e+00 4.60e-07R  5
16r 2.5464204e-02 3.56e+01 9.99e+02  -0.4 5.23e+02    -  1.19e-03 5.15e-04f  1
17  2.5464199e-02 3.56e+01 1.84e+01  -1.0 1.70e+02    -  7.36e-05 1.10e-05h  1
18  2.5464182e-02 3.56e+01 5.50e+01  -1.0 5.87e+02    -  7.29e-06 3.93e-05f  1
19r 2.5464182e-02 3.56e+01 9.99e+02  -0.6 0.00e+00    -  0.00e+00 4.49e-07R  4
iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
20r 2.7129142e-02 3.30e+01 9.98e+02  -0.6 2.08e+02    -  1.32e-03 8.68e-05f  1
21r 5.1678170e-02 1.20e+00 9.95e+02  -0.6 9.59e+01    -  3.38e-03 2.31e-03f  1
22r 5.1678170e-02 1.20e+00 9.99e+02  -1.0 0.00e+00    -  0.00e+00 3.15e-07R  3
23r 5.6729089e-02 1.18e+00 9.98e+02  -1.0 3.56e+01    -  1.77e-03 6.40e-04f  1
24r 6.6002832e-02 1.14e+00 9.97e+02  -1.0 1.42e+01    -  1.58e-03 1.20e-03f  1
25r 9.8505667e-02 9.98e-01 9.92e+02  -1.0 7.28e+00    -  5.20e-03 4.47e-03f  1
26r 1.0225287e-01 9.87e-01 9.92e+02  -1.0 4.78e+00    -  5.96e-03 7.89e-04f  1
27r 1.1636717e-01 9.44e-01 9.88e+02  -1.0 4.45e+00    -  3.90e-03 3.22e-03f  1
28r 1.3177033e-01 8.84e-01 9.84e+02  -1.0 3.30e+00    -  4.43e-03 4.73e-03f  1
29r 1.4160321e-01 8.19e-01 9.79e+02  -1.0 2.44e+00    -  5.28e-03 5.24e-03f  1
iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
30  1.4145724e-01 8.19e-01 2.55e+04  -1.0 5.48e+02    -  6.79e-05 2.14e-04f  1
31  1.4096248e-01 8.19e-01 2.34e+04  -1.0 4.57e+02    -  1.35e-04 8.43e-04f  1
32  1.4095301e-01 8.19e-01 2.34e+04  -1.0 3.28e+02    -  1.40e-04 1.17e-05h  1
33  1.4078924e-01 8.19e-01 3.61e+05  -1.0 1.79e+03    -  3.25e-06 1.77e-04f  1
34  1.4056815e-01 8.19e-01 5.35e+05  -1.0 9.60e+02    -  7.83e-05 2.49e-04f  1
35  1.4037264e-01 8.19e-01 5.58e+05  -1.0 1.30e+03    -  7.61e-05 2.09e-04h  1
36  1.4023681e-01 8.19e-01 7.31e+05  -1.0 2.47e+03    -  4.95e-05 1.32e-04h  1
37  1.4016510e-01 8.19e-01 7.48e+05  -1.0 3.86e+03    -  3.59e-05 5.03e-05h  1
38  1.4009700e-01 8.19e-01 1.04e+06  -1.0 6.74e+03    -  1.70e-05 3.96e-05h  1
39  1.4005369e-01 8.19e-01 1.20e+06  -1.0 1.58e+04    -  7.56e-06 1.49e-05h  1
iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
40r 1.4005369e-01 8.19e-01 1.00e+03  -1.0 0.00e+00    -  0.00e+00 2.96e-07R  6
41r 1.4136641e-01 7.99e-01 4.09e+03  -1.0 1.76e+00    -  9.82e-03 1.98e-03f  1
42r 1.4730479e-01 6.38e-01 6.95e+03  -1.0 3.10e+00    -  3.85e-03 1.35e-02f  1
43r 1.4730479e-01 6.38e-01 9.98e+02  -1.0 0.00e+00    -  0.00e+00 1.15e-07R  2
44r 1.4779590e-01 5.25e-01 9.92e+02  -1.0 1.79e+00    -  1.24e-02 6.28e-03f  1
45r 1.4843078e-01 3.49e-01 2.60e+03  -1.0 2.85e+00    -  6.70e-03 2.53e-02f  1
46r 1.4782353e-01 2.71e-01 1.45e+03  -1.0 2.66e+00    -  8.47e-03 1.73e-02f  1
47r 1.4758686e-01 2.52e-01 1.50e+03  -1.0 2.80e+00    -  8.44e-03 6.06e-03f  1
48r 1.4745003e-01 2.44e-01 1.51e+03  -1.0 3.63e+00    -  6.47e-03 3.05e-03f  1
49r 1.4717730e-01 2.33e-01 1.56e+03  -1.0 4.63e+00    -  5.01e-03 4.92e-03f  1
iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
50r 1.4678673e-01 2.20e-01 3.74e+03  -1.0 4.94e+00    -  5.58e-03 6.43e-03f  1
51r 1.4575626e-01 1.88e-01 1.47e+04  -1.0 4.72e+00    -  6.02e-03 1.71e-02f  1
52r 1.4566201e-01 1.80e-01 1.47e+04  -1.0 2.25e+00    -  6.27e-03 2.07e-03f  1
53r 1.4548158e-01 1.70e-01 3.55e+04  -1.0 1.39e+00   2.0 5.53e-03 6.91e-03f  1
54r 1.4539014e-01 1.27e-01 2.97e+04  -1.0 2.34e-02   4.2 1.06e-01 3.15e-02f  1
55r 1.4533095e-01 1.18e-01 2.64e+04  -1.0 1.08e-01   3.8 2.99e-02 1.27e-02f  1
56r 1.4512869e-01 7.71e-02 2.27e+04  -1.0 1.78e-02   4.2 1.28e-01 7.20e-02f  1
57  1.4490266e-01 7.71e-02 7.08e+03  -1.0 1.54e+01    -  5.20e-03 5.83e-04f  1
58  1.4452813e-01 7.71e-02 2.91e+05  -1.0 6.73e+02    -  1.48e-05 4.33e-04h  1
59  1.4422222e-01 7.71e-02 4.49e+05  -1.0 6.49e+02    -  2.53e-04 2.27e-04h  1
iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
60  1.4370807e-01 7.71e-02 5.83e+05  -1.0 1.27e+03    -  1.31e-04 2.21e-04h  1
61  1.4297321e-01 7.71e-02 5.92e+05  -1.0 7.02e+02    -  6.73e-05 1.25e-04h  1
62  1.4072428e-01 7.71e-02 1.36e+06  -1.0 1.06e+03    -  1.31e-04 2.95e-04h  1
63  1.2873821e-01 7.71e-02 4.49e+06  -1.0 1.04e+02    -  4.09e-04 1.68e-03f  1
64  1.2871221e-01 7.71e-02 4.49e+06  -1.0 1.21e+02    -  4.74e-03 3.50e-05h  1
65  1.2490836e-01 7.70e-02 7.36e+06  -1.0 4.20e+01    -  3.42e-03 4.95e-03f  1
66  8.4787480e-02 7.64e-02 1.64e+07  -1.0 5.03e+00    -  4.98e-03 8.35e-02f  1
67  7.8553840e-02 7.59e-02 1.47e+07  -1.0 6.66e+00    -  6.66e-02 6.69e-02f  1
68  7.5543485e-02 7.56e-02 1.40e+07  -1.0 6.37e+00    -  1.49e-01 3.70e-02f  1
69  7.5515373e-02 7.56e-02 3.49e+08  -1.0 6.19e+00    -  1.45e-01 3.64e-04h  1
iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
70  7.5512656e-02 7.56e-02 3.79e+11  -1.0 1.13e+01    -  3.81e-02 3.53e-05h  1
71  7.5504268e-02 7.56e-02 6.40e+13  -1.0 1.57e+01    -  1.86e-02 1.09e-04h  1
72r 7.5504268e-02 7.56e-02 1.00e+03  -1.0 0.00e+00    -  0.00e+00 3.41e-07R  3
73r 7.5563622e-02 7.53e-02 9.89e+03  -1.0 1.63e+01    -  1.05e-02 1.64e-03f  1
74r 7.5792890e-02 7.31e-02 8.75e+03  -1.0 1.43e+01    -  5.81e-03 6.45e-03f  1
75  7.3203256e-02 7.29e-02 4.82e+05  -1.0 6.13e+00    -  1.72e-01 3.35e-02f  1
76  7.3162034e-02 7.29e-02 2.62e+08  -1.0 5.81e+00    -  1.83e-01 5.56e-04h  1
77  7.2137141e-02 7.27e-02 5.87e+09  -1.0 5.61e+00    -  2.15e-01 1.39e-02h  1
78  7.2100734e-02 7.27e-02 1.74e+12  -1.0 5.52e+00    -  1.40e-01 5.00e-04h  1
79  7.2099621e-02 7.27e-02 2.64e+16  -1.0 5.52e+00    -  2.34e-01 1.53e-05h  1
iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
80r 7.2099621e-02 7.27e-02 1.00e+03  -1.0 0.00e+00    -  0.00e+00 3.98e-07R  3
81r 7.2182761e-02 7.15e-02 5.54e+03  -1.0 1.77e+01    -  7.29e-03 3.50e-03f  1
82  7.1738248e-02 7.15e-02 1.79e+07  -1.0 5.74e+00    -  2.21e-01 6.10e-03f  1
83  6.5183013e-02 7.07e-02 3.80e+08  -1.0 5.68e+00    -  3.96e-03 9.06e-02f  1
84  6.5116414e-02 7.07e-02 3.91e+08  -1.0 5.20e+00    -  9.29e-02 1.02e-03h  1
85  6.5115747e-02 7.07e-02 3.79e+08  -1.0 5.18e+00    -  1.41e-02 1.02e-05h  1
86r 6.5115747e-02 7.07e-02 1.00e+03  -1.0 0.00e+00    -  0.00e+00 4.93e-07R  4
87r 6.5221569e-02 6.95e-02 6.98e+03  -1.0 2.05e+01    -  8.26e-03 2.34e-03f  1
88  6.5182297e-02 6.95e-02 1.16e+07  -1.0 4.99e+00    -  3.10e-01 6.03e-04h  1
89r 6.5182297e-02 6.95e-02 1.00e+03  -1.0 0.00e+00    -  0.00e+00 3.78e-07R  5
iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
90r 6.5319853e-02 6.80e-02 5.63e+03  -1.0 1.96e+01    -  8.66e-03 3.07e-03f  1
91  6.5266535e-02 6.80e-02 1.37e+07  -1.0 4.90e+00    -  3.20e-01 8.18e-04h  1
92r 6.5266535e-02 6.80e-02 1.00e+03  -1.0 0.00e+00    -  0.00e+00 2.56e-07R  6
93r 6.5440982e-02 6.57e-02 4.73e+03  -1.0 1.85e+01    -  9.11e-03 4.01e-03f  1
94  6.5368738e-02 6.57e-02 9.66e+06  -1.0 4.81e+00    -  2.94e-01 1.11e-03h  1
95  6.5368015e-02 6.57e-02 2.36e+09  -1.0 4.80e+00    -  3.61e-01 1.11e-05h  1
96r 6.5368015e-02 6.57e-02 1.00e+03  -1.0 0.00e+00    -  0.00e+00 3.08e-07R  5
97r 6.5504081e-02 6.31e-02 4.46e+03  -1.0 1.73e+01    -  1.05e-02 4.58e-03f  1
98  6.5418433e-02 6.31e-02 6.84e+06  -1.0 4.71e+00    -  2.72e-01 1.31e-03h  1
99  6.5417575e-02 6.31e-02 1.64e+09  -1.0 4.76e+00    -  3.63e-01 1.32e-05h  1
iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
100r 6.5417575e-02 6.31e-02 1.00e+03  -1.0 0.00e+00    -  0.00e+00 3.48e-07R  5
101r 6.5619873e-02 5.94e-02 4.52e+03  -1.0 1.58e+01    -  1.35e-02 6.84e-03f  1
102r 6.6453038e-02 4.74e-02 1.11e+04  -1.0 1.17e+01    -  4.83e-02 2.89e-02f  1
103r 6.7406490e-02 7.23e-02 9.49e+03  -1.0 4.55e+00    -  2.76e-02 3.52e-02f  1
104r 6.7340988e-02 7.22e-02 2.76e+04  -1.0 3.44e-02   4.0 8.50e-02 3.51e-02f  1
105r 6.7352639e-02 7.17e-02 2.73e+04  -1.0 4.13e+00    -  8.44e-04 4.11e-04f  1
106r 6.7213982e-02 7.15e-02 1.56e+04  -1.0 6.13e-02   3.5 2.79e-02 5.72e-02f  1
107r 6.7187230e-02 7.15e-02 2.61e+04  -1.0 5.92e-02   3.9 4.76e-02 1.52e-02f  1
108r 6.7106847e-02 7.14e-02 3.45e+04  -1.0 4.31e-02   3.5 6.87e-02 3.47e-02f  1
109r 6.7066182e-02 7.13e-02 3.33e+04  -1.0 1.50e-01   3.9 2.16e-02 2.35e-02f  1
iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
110r 6.7006506e-02 7.13e-02 3.24e+04  -1.0 4.01e-02   3.4 2.69e-02 2.68e-02f  1
111r 6.6859346e-02 7.12e-02 6.05e+04  -1.0 1.29e-02   3.8 2.01e-01 8.60e-02f  1
112r 6.6575748e-02 7.09e-02 4.81e+04  -1.0 8.58e-03   4.3 2.85e-01 2.68e-01f  1
113r 6.6559828e-02 7.09e-02 4.82e+04  -1.0 1.72e-01   3.8 2.11e-02 1.16e-02f  1
114r 6.6364881e-02 7.06e-02 7.22e+04  -1.0 8.35e-03   4.2 5.17e-01 2.14e-01f  1
115r 6.6121763e-02 7.02e-02 6.05e+04  -1.0 1.73e-03   4.6 1.00e+00 5.22e-01f  1
116r 6.5586391e-02 6.86e-02 1.03e+04  -1.0 2.55e-03   4.2 1.00e+00 8.02e-01f  1
117r 6.4949064e-02 6.44e-02 9.59e+03  -1.0 4.71e-03   3.7 1.00e+00 1.00e+00f  1
118r 6.4885306e-02 5.76e-02 3.22e+04  -1.0 1.47e-02   3.2 1.00e+00 7.92e-01f  1
119r 6.4897911e-02 5.40e-02 7.87e+03  -1.0 3.55e-03   3.6 1.00e+00 1.00e+00f  1
iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
120r 6.5282534e-02 4.66e-02 1.16e+04  -1.0 7.48e-03   3.2 1.00e+00 1.00e+00f  1
121r 6.7015934e-02 3.28e-02 7.20e+03  -1.0 1.37e-02   2.7 1.00e+00 1.00e+00f  1
122r 7.2635183e-02 2.45e-02 1.05e+04  -1.0 2.59e-02   2.2 1.00e+00 1.00e+00f  1
123r 7.3336436e-02 3.01e-02 3.51e+04  -1.0 9.95e-01   1.7 5.47e-02 6.05e-02F  1
124r 7.4118433e-02 3.05e-02 1.88e+04  -1.0 5.77e-03   3.1 1.00e+00 1.00e+00f  1
125r 7.5881054e-02 2.96e-02 1.00e+04  -1.0 1.32e-02   2.6 9.48e-01 7.60e-01f  1
126r 7.5952551e-02 2.96e-02 2.88e+03  -1.0 1.86e-03   3.9 1.00e+00 6.56e-01f  1
127r 7.6279414e-02 2.94e-02 8.69e+02  -1.0 1.85e-03   3.4 1.00e+00 1.00e+00f  1
128r 7.7252002e-02 2.88e-02 5.05e+02  -1.0 5.53e-03   3.0 1.00e+00 1.00e+00f  1
129r 8.0052578e-02 2.79e-02 2.40e+03  -1.0 1.67e-02   2.5 1.00e+00 1.00e+00f  1
iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
130r 8.2021764e-02 2.83e-02 8.20e+04  -1.0 1.13e-01   2.0 3.19e-01 2.77e-01F  1
131r 8.4876875e-02 2.83e-02 1.58e+04  -1.0 1.90e-02   2.4 1.00e+00 1.00e+00f  1
132  8.0249864e-02 2.67e-02 1.25e+04  -1.0 2.50e+00    -  7.13e-02 5.54e-02h  1
133  7.4636193e-02 2.45e-02 3.46e+04  -1.0 3.00e+00    -  8.77e-02 7.12e-02h  1
134  7.3809769e-02 2.41e-02 7.25e+05  -1.0 2.49e+00    -  1.74e-01 1.13e-02h  1
135  7.3798438e-02 2.41e-02 5.19e+07  -1.0 7.94e+00    -  9.89e-03 1.56e-04h  1
136  7.3620037e-02 2.41e-02 1.10e+07  -1.0 1.49e+02    -  6.02e-04 2.45e-03h  1
137  7.3605855e-02 2.40e-02 2.11e+09  -1.0 9.24e+01    -  4.43e-03 1.95e-04h  1
138  7.3599343e-02 2.40e-02 2.14e+10  -1.0 4.53e+02    -  8.43e-04 8.97e-05h  1
139r 7.3599343e-02 2.40e-02 1.00e+03  -1.0 0.00e+00    -  0.00e+00 4.49e-07R  4
iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
140r 7.4408520e-02 2.43e-02 1.01e+04  -1.0 5.67e-01    -  6.63e-02 1.97e-02f  1
141r 7.7821300e-02 3.19e-02 1.31e+04  -1.0 4.34e-01    -  1.27e-01 8.13e-02f  1
142  7.7333980e-02 3.17e-02 5.09e+05  -1.0 2.07e+00    -  1.28e-01 6.37e-03f  1
143  7.7291431e-02 3.17e-02 1.49e+08  -1.0 3.55e+00    -  1.10e-01 5.59e-04h  1
144  7.6515671e-02 3.14e-02 2.88e+09  -1.0 3.36e+00    -  1.82e-01 1.02e-02h  1
145  7.6487635e-02 3.14e-02 1.20e+11  -1.0 2.92e+01    -  1.48e-02 3.73e-04h  1
146r 7.6487635e-02 3.14e-02 1.00e+03  -1.0 0.00e+00    -  0.00e+00 4.25e-07R  5
147r 7.8138416e-02 3.54e-02 1.80e+04  -1.0 4.35e-01    -  1.71e-01 4.45e-02f  1
148r 8.3135749e-02 4.69e-02 7.66e+04  -1.0 4.92e-01    -  7.66e-01 1.27e-01f  1
149  7.9891006e-02 4.59e-02 2.42e+05  -1.0 3.41e+00    -  8.57e-02 3.96e-02f  1
iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
150  7.8120496e-02 4.53e-02 1.72e+07  -1.0 2.43e+00    -  5.26e-01 2.25e-02h  1
151  7.8075810e-02 4.53e-02 7.19e+09  -1.0 2.35e+00    -  2.22e-01 5.81e-04h  1
152  7.7942727e-02 4.52e-02 8.56e+11  -1.0 2.28e+00    -  2.06e-01 1.73e-03h  1
153  7.7940793e-02 4.52e-02 2.89e+16  -1.0 2.72e+00    -  8.52e-01 2.52e-05h  1
154r 7.7940793e-02 4.52e-02 1.00e+03  -1.0 0.00e+00    -  0.00e+00 1.40e-07R  2
155r 7.8114585e-02 4.55e-02 6.99e+04  -1.0 1.59e+00    -  2.28e-01 9.50e-03f  3
156r 8.3677416e-02 6.22e-02 8.70e+04  -1.0 1.58e+00    -  9.92e-01 2.53e-01f  1
157  8.3635098e-02 6.22e-02 2.59e+06  -1.0 2.04e+00    -  5.89e-01 5.12e-04f  5
158  8.2999202e-02 6.20e-02 2.51e+08  -1.0 2.09e+00    -  7.47e-01 7.70e-03h  1
159  8.1933530e-02 6.16e-02 6.76e+09  -1.0 1.18e+00    -  3.47e-01 1.30e-02h  1
iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
160  8.1807679e-02 6.16e-02 3.37e+12  -1.0 1.66e+00    -  7.56e-01 1.56e-03h  1
161  8.1805967e-02 6.16e-02 2.54e+16  -1.0 2.93e+00    -  1.61e-01 2.12e-05h  1
162  8.1805967e-02 6.16e-02 2.51e+14  -1.0 6.77e-09  16.0 9.90e-01 1.00e+00h  1
163r 8.1805967e-02 6.16e-02 1.00e+03  -1.0 0.00e+00    -  0.00e+00 1.17e-07R  2
164r 8.6400165e-02 6.59e-02 1.10e+05  -1.0 1.39e+00    -  6.83e-01 1.21e-01f  1
165r 1.2303118e-01 8.53e-02 3.10e+04  -1.0 1.16e+00    -  8.31e-01 7.64e-01f  1
166r 1.2359049e-01 9.10e-02 3.79e+04  -1.0 1.36e-02   2.0 8.58e-01 3.33e-01f  1
167r 1.2396345e-01 8.75e-02 1.06e+04  -1.0 1.38e-02   2.4 1.00e+00 5.88e-01f  1
168r 1.2575264e-01 7.05e-02 4.85e+02  -1.0 2.22e-02   1.9 1.00e+00 1.00e+00f  1
169r 1.2983673e-01 3.79e-02 3.03e+02  -1.0 5.61e-02   1.5 1.00e+00 1.00e+00f  1
iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
170r 1.2989395e-01 3.80e-02 3.13e+03  -1.0 3.45e+00   1.0 2.91e-02 8.66e-03f  3
171r 1.2996510e-01 3.76e-02 6.80e+02  -1.0 1.11e-03   3.2 1.00e+00 1.00e+00f  1
172r 1.3017631e-01 3.71e-02 7.92e+01  -1.0 2.29e-03   2.7 1.00e+00 1.00e+00f  1
173r 1.3078032e-01 3.29e-02 1.76e+01  -1.0 6.28e-03   2.3 1.00e+00 1.00e+00f  1
174r 1.3234891e-01 2.34e-02 5.98e+02  -1.0 3.11e-02   1.8 1.00e+00 1.00e+00f  1
175r 1.3259492e-01 2.27e-02 3.69e+04  -1.0 4.00e-01   1.3 3.10e-01 7.58e-02f  2
176r 1.3406520e-01 1.40e-02 3.84e+01  -1.0 3.73e-02   1.7 1.00e+00 1.00e+00f  1
177r 1.3697540e-01 1.21e-02 1.90e+02  -1.0 8.43e-02   1.3 1.00e+00 1.00e+00f  1
178r 1.3790207e-01 1.20e-02 6.40e+01  -1.0 1.42e-02   1.7 1.00e+00 1.00e+00f  1
179r 1.3969952e-01 1.32e-02 8.38e+02  -1.0 4.38e-02   1.2 1.00e+00 1.00e+00f  1
iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
180r 1.4026468e-01 1.39e-02 2.09e+02  -1.0 1.94e-02   1.6 1.00e+00 1.00e+00f  1
181r 1.4046219e-01 1.43e-02 4.53e+01  -1.0 7.58e-03   2.1 1.00e+00 1.00e+00f  1
182r 1.4095097e-01 1.56e-02 1.16e+03  -1.0 4.18e-02   1.6 1.00e+00 1.00e+00f  1
183r 1.4112123e-01 1.63e-02 6.22e+01  -1.0 2.16e-02   2.0 1.00e+00 1.00e+00f  1
184r 1.4153884e-01 1.81e-02 3.22e+02  -1.0 6.21e-02   1.5 1.00e+00 1.00e+00f  1
185r 1.4168358e-01 1.88e-02 6.74e+01  -1.0 1.92e-02   2.0 1.00e+00 1.00e+00f  1
186r 1.4183622e-01 2.03e-02 3.63e+03  -1.0 1.00e-01   1.5 6.33e-01 4.30e-01f  2
187r 1.4197415e-01 2.12e-02 3.52e+01  -1.0 2.48e-02   1.9 1.00e+00 1.00e+00f  1
188r 1.4229977e-01 2.50e-02 2.57e+02  -1.0 6.82e-02   1.4 1.00e+00 1.00e+00f  1
189r 1.4263582e-01 5.60e-02 4.48e+03  -1.0 1.04e-01   1.0 1.00e+00 7.62e-01F  1
iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
190r 1.4265834e-01 6.84e-02 1.68e+03  -1.0 1.53e-02   2.3 1.00e+00 8.85e-01h  1
191r 1.4272531e-01 7.28e-02 8.14e+01  -1.0 2.52e-02   1.8 1.00e+00 1.00e+00f  1
192r 1.4283214e-01 1.08e-01 1.35e+03  -1.0 6.73e-02   1.3 1.00e+00 1.00e+00f  1
193r 1.4281762e-01 1.07e-01 9.84e+03  -1.0 8.82e-01   0.9 1.29e-01 4.01e-02f  3
194r 1.4280315e-01 1.19e-01 1.36e+03  -1.0 4.68e-02   1.3 1.00e+00 1.00e+00f  1
195r 1.4254558e-01 1.11e-01 1.11e+03  -1.0 8.02e-02   0.8 1.00e+00 1.00e+00f  1
196r 1.4244677e-01 1.08e-01 4.27e+01  -1.0 1.80e-02   1.2 1.00e+00 1.00e+00f  1
197r 1.4226538e-01 1.04e-01 1.32e+02  -1.0 4.53e-02   0.8 1.00e+00 1.00e+00f  1
198r 1.4221142e-01 1.01e-01 5.69e+01  -1.0 1.43e-02   1.2 1.00e+00 1.00e+00f  1
199r 1.4220172e-01 1.01e-01 4.59e+03  -1.0 9.52e-02   0.7 1.00e+00 2.50e-01f  3
iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
200r 1.4220122e-01 1.00e-01 1.71e+04  -1.0 8.56e-02   1.1 1.00e+00 1.25e-01f  4
201r 1.4218741e-01 9.95e-02 2.91e+00  -1.0 7.11e-03   1.6 1.00e+00 1.00e+00f  1
202r 1.2582161e-01 1.03e-01 1.30e+04  -1.7 1.46e-01   1.1 3.74e-01 3.62e-01f  1
203r 1.1051968e-01 1.15e-01 3.68e+04  -1.7 3.54e-02   1.5 8.71e-01 6.59e-01f  1
204r 1.0129625e-01 1.44e-01 2.96e+04  -1.7 1.57e-02   1.9 1.00e+00 8.47e-01f  1
205r 9.7073006e-02 1.42e-01 2.51e+04  -1.7 1.73e-01   1.5 2.22e-01 1.83e-01f  1
206r 8.7124928e-02 1.35e-01 8.58e+03  -1.7 1.16e-02   1.9 1.00e+00 1.00e+00f  1
207r 8.3201403e-02 1.23e-01 7.37e+03  -1.7 1.20e-01   1.4 2.49e-01 2.11e-01f  1
208r 7.5093725e-02 8.20e-02 2.73e+03  -1.7 1.63e-02   1.8 1.00e+00 1.00e+00f  1
209r 6.9950640e-02 3.25e-02 2.90e+03  -1.7 8.71e-02   1.4 4.35e-01 3.41e-01f  1
iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
210r 6.4656153e-02 7.45e-03 2.10e+03  -1.7 1.43e-02   1.8 1.00e+00 7.91e-01f  1
246  2.8072332e-03 1.66e-05 2.49e+06  -1.0 2.22e-01    -  1.00e+00 4.81e-02h  5
247  2.8349267e-03 1.57e-05 2.70e+06  -1.0 2.04e-01    -  1.00e+00 4.84e-02h  5
248  2.8629932e-03 1.50e-05 2.90e+06  -1.0 1.97e-01    -  1.00e+00 4.84e-02h  5
249  2.8917742e-03 1.43e-05 3.07e+06  -1.0 1.89e-01    -  1.00e+00 4.84e-02h  5
iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
250  3.3664069e-03 6.69e-05 1.57e+06  -1.0 1.81e-01    -  1.00e+00 7.75e-01w  1
251  5.5300533e-03 1.11e-03 1.02e+07  -1.0 1.30e+00    -  2.62e-01 9.87e-01w  1
252  1.0629220e-02 3.18e-03 3.01e+06  -1.0 1.57e+00    -  3.65e-01 1.00e+00w  1
253  2.9511033e-03 1.29e-05 3.06e+06  -1.0 3.26e-01    -  1.00e+00 9.68e-02h  3
254  2.9848641e-03 1.23e-05 3.84e+06  -1.0 1.66e-01    -  1.00e+00 4.81e-02h  5
255  3.0190653e-03 1.17e-05 4.72e+06  -1.0 1.54e-01    -  1.00e+00 4.87e-02h  5
256  3.0542668e-03 1.11e-05 5.78e+06  -1.0 1.47e-01    -  1.00e+00 4.89e-02h  5
257  3.0906899e-03 1.05e-05 7.03e+06  -1.0 1.39e-01    -  1.00e+00 4.91e-02h  5
258  3.1095384e-03 1.03e-05 8.72e+06  -1.0 1.32e-01    -  1.00e+00 2.47e-02h  6
259  3.1468402e-03 9.78e-06 1.04e+07  -1.0 1.29e-01    -  1.00e+00 4.97e-02h  5
iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
260  3.1665035e-03 9.53e-06 1.28e+07  -1.0 1.38e-01    -  1.00e+00 2.49e-02h  6
261  3.2055680e-03 9.07e-06 1.51e+07  -1.0 1.24e-01    -  1.00e+00 5.02e-02h  5
262  3.2262213e-03 8.85e-06 1.83e+07  -1.0 1.60e-01    -  1.00e+00 2.52e-02h  6
263  3.8835449e-03 4.04e-05 7.93e+06  -1.0 1.48e-01    -  1.00e+00 8.14e-01w  1
264  5.8638686e-03 1.20e-03 1.62e+08  -1.0 1.17e+00    -  1.19e-01 1.00e+00w  1
265  9.7329055e-03 2.80e-03 2.36e+06  -1.0 1.35e+00    -  8.82e-01 1.00e+00w  1
266  3.2467627e-03 8.62e-06 2.19e+07  -1.0 1.50e+00    -  1.00e+00 2.54e-02h  5
267  3.2675339e-03 8.41e-06 2.60e+07  -1.0 1.51e-01    -  1.00e+00 2.55e-02h  6
268  3.2886208e-03 8.20e-06 3.08e+07  -1.0 1.59e-01    -  1.00e+00 2.57e-02h  6
269  3.3100455e-03 7.99e-06 3.63e+07  -1.0 1.67e-01    -  1.00e+00 2.58e-02h  6
iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
270  3.3318361e-03 7.78e-06 4.25e+07  -1.0 1.77e-01    -  1.00e+00 2.59e-02h  6
312  4.0689751e-03 4.78e-06 2.09e+08  -1.0 6.42e-01    -  1.00e+00 1.56e-02h  7
313  4.0906299e-03 4.81e-06 1.95e+08  -1.0 6.55e-01    -  1.00e+00 1.56e-02h  7
314  4.1126371e-03 4.83e-06 1.80e+08  -1.0 6.68e-01    -  1.00e+00 1.56e-02h  7
315  5.5446739e-03 4.27e-04 3.22e+07  -1.0 6.81e-01    -  1.00e+00 1.00e+00w  1
316  1.0125712e-02 4.28e-03 1.16e+09  -1.0 1.95e+00    -  1.95e-01 1.00e+00w  1
317  1.6101826e-02 1.38e-04 1.31e+06  -1.0 2.15e-01    -  1.00e+00 1.00e+00w  1
318  4.1350127e-03 4.86e-06 1.64e+08  -1.0 8.33e-01    -  1.00e+00 1.56e-02h  6
319  4.1577748e-03 4.89e-06 1.48e+08  -1.0 6.94e-01    -  1.00e+00 1.56e-02h  7
iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
320  4.1693598e-03 4.88e-06 1.33e+08  -1.0 7.08e-01    -  1.00e+00 7.81e-03h  8
321  4.1810345e-03 4.87e-06 1.17e+08  -1.0 7.12e-01    -  1.00e+00 7.81e-03h  8
322  4.1928460e-03 4.86e-06 1.03e+08  -1.0 7.20e-01    -  1.00e+00 7.81e-03h  8
323  4.2048136e-03 4.85e-06 8.88e+07  -1.0 7.28e-01    -  1.00e+00 7.81e-03h  8
324  4.2169515e-03 4.85e-06 7.57e+07  -1.0 7.37e-01    -  1.00e+00 7.81e-03h  8
325  4.2292851e-03 4.84e-06 6.35e+07  -1.0 7.47e-01    -  1.00e+00 7.81e-03h  8
326  4.2418574e-03 4.83e-06 5.25e+07  -1.0 7.59e-01    -  1.00e+00 7.81e-03h  8
327  4.2547649e-03 4.82e-06 4.32e+07  -1.0 7.73e-01    -  1.00e+00 7.81e-03h  8
328  5.9950385e-03 5.45e-04 5.75e+06  -1.0 8.00e-01    -  7.71e-01 1.00e+00w  1
329  6.6127800e-03 5.45e-04 1.76e+07  -1.0 8.41e+01    -  3.70e-04 4.93e-03w  1
iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
330  7.5349653e-03 4.38e-07 1.85e+17  -1.0 5.32e-03  15.5 1.53e-01 1.00e+00h  1
MUMPS returned INFO(1) = -9 and requires more memory, reallocating.  Attempt 1
Increasing icntl[13] from 1000 to 2000.
MUMPS returned INFO(1) = -9 and requires more memory, reallocating.  Attempt 2
Increasing icntl[13] from 2000 to 4000.
331  7.5575806e-03 2.73e-09 8.70e+15  -1.0 7.81e-04  15.0 1.00e+00 1.00e+00f  1
332  7.5575628e-03 7.99e-15 3.78e+11  -1.0 1.28e-06  14.6 1.00e+00 1.00e+00h  1
333  7.5575628e-03 4.44e-16 2.79e+02  -1.7 2.26e-12  14.1 1.00e+00 1.00e+00h  1
334  7.5575628e-03 4.44e-16 4.85e+00  -5.7 1.18e-13  13.6 1.00e+00 1.00e+00h  1
335  7.5575628e-03 4.44e-16 2.50e-02  -5.7 1.83e-15  13.1 1.00e+00 1.00e+00   0
336  7.5575628e-03 4.44e-16 4.61e-03  -5.7 1.01e-15  12.7 1.00e+00 1.00e+00   0
337  7.5575628e-03 4.44e-16 5.70e-04  -5.7 3.74e-16  12.2 1.00e+00 1.00e+00   0
338  7.5575628e-03 4.44e-16 5.44e-04  -5.7 1.07e-15  11.7 1.00e+00 1.00e+00   0
339  7.5575628e-03 4.44e-16 5.17e-04  -5.7 3.05e-15  11.2 1.00e+00 1.00e+00   0
iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
340  7.5575628e-03 4.44e-16 5.25e-04  -5.7 9.31e-15  10.8 1.00e+00 1.00e+00h  1
341  7.5575628e-03 4.44e-16 5.29e-04  -5.7 2.81e-14  10.3 1.00e+00 1.00e+00h  1
342  7.5575628e-03 4.44e-16 5.26e-04  -5.7 8.38e-14   9.8 1.00e+00 1.00e+00h  1
343  7.5575628e-03 4.44e-16 5.24e-04  -5.7 2.51e-13   9.3 1.00e+00 1.00e+00h  1
344  7.5575628e-03 4.44e-16 5.25e-04  -5.7 7.53e-13   8.8 1.00e+00 1.00e+00h  1
345  7.5575628e-03 4.44e-16 5.25e-04  -5.7 2.26e-12   8.4 1.00e+00 1.00e+00h  1
346  7.5575628e-03 4.44e-16 5.25e-04  -5.7 6.78e-12   7.9 1.00e+00 1.00e+00h  1
347  7.5575628e-03 4.44e-16 5.25e-04  -5.7 2.03e-11   7.4 1.00e+00 1.00e+00h  1
348  7.5575628e-03 4.44e-16 5.25e-04  -5.7 6.10e-11   6.9 1.00e+00 1.00e+00h  1
349  7.5575628e-03 4.44e-16 5.25e-04  -5.7 1.83e-10   6.5 1.00e+00 1.00e+00h  1
iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
350  7.5575628e-03 4.44e-16 5.25e-04  -5.7 5.49e-10   6.0 1.00e+00 1.00e+00h  1
351  7.5575628e-03 4.44e-16 5.25e-04  -5.7 1.65e-09   5.5 1.00e+00 1.00e+00h  1
352  7.5575628e-03 4.44e-16 5.25e-04  -5.7 4.94e-09   5.0 1.00e+00 1.00e+00h  1
353  7.5575628e-03 4.44e-16 5.25e-04  -5.7 1.48e-08   4.5 1.00e+00 1.00e+00h  1
354  7.5575629e-03 4.44e-16 5.25e-04  -5.7 4.45e-08   4.1 1.00e+00 1.00e+00h  1
355  7.5575629e-03 4.44e-16 5.25e-04  -5.7 1.33e-07   3.6 1.00e+00 1.00e+00h  1
356  7.5575630e-03 4.44e-16 5.25e-04  -5.7 4.00e-07   3.1 1.00e+00 1.00e+00h  1
357  7.5575633e-03 5.69e-16 5.25e-04  -5.7 1.20e-06   2.6 1.00e+00 1.00e+00h  1
358  7.5575641e-03 5.18e-15 5.24e-04  -5.7 3.60e-06   2.2 1.00e+00 1.00e+00h  1
359  7.5575667e-03 4.64e-14 5.23e-04  -5.7 1.08e-05   1.7 1.00e+00 1.00e+00h  1
iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
360  7.5575746e-03 4.15e-13 5.21e-04  -5.7 3.22e-05   1.2 1.00e+00 1.00e+00h  1
361  7.5575981e-03 3.67e-12 5.12e-04  -5.7 9.49e-05   0.7 1.00e+00 1.00e+00h  1
362  7.5576693e-03 3.14e-11 4.89e-04  -5.7 2.72e-04   0.3 1.00e+00 1.00e+00h  1
363  7.5578876e-03 2.47e-10 4.34e-04  -5.7 7.24e-04  -0.2 1.00e+00 1.00e+00h  1
364  7.5585735e-03 1.68e-09 1.84e-03  -5.7 1.70e-03  -0.7 1.00e+00 1.00e+00h  1
365  7.5607791e-03 9.50e-09 6.74e-03  -5.7 3.56e-03  -1.2 1.00e+00 1.00e+00h  1
366  7.5678705e-03 4.57e-08 1.93e-02  -5.7 7.12e-03  -1.7 1.00e+00 1.00e+00h  1
367  7.5902322e-03 2.24e-07 4.07e-02  -5.7 1.39e-02  -2.1 1.00e+00 1.00e+00h  1
368  7.6594408e-03 2.05e-06 7.09e-02  -5.7 2.66e-02  -2.6 1.00e+00 1.00e+00h  1
369  7.8736336e-03 1.94e-05 1.17e-01  -5.7 7.37e-02  -3.1 1.00e+00 1.00e+00h  1
iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
370  8.5895324e-03 2.15e-04 1.94e-01  -5.7 2.39e-01  -3.6 1.00e+00 1.00e+00h  1
371  1.0568509e-02 1.63e-03 2.92e+00  -5.7 1.13e+00  -4.0 9.01e-01 5.25e-01h  1
372  1.3005687e-02 7.79e-04 4.97e-01  -5.7 3.94e-01  -3.6 1.00e+00 1.00e+00f  1
373  1.6306827e-02 1.60e-03 6.32e+00  -5.7 3.32e+00  -4.1 3.45e-01 1.38e-01h  1
374  2.0759846e-02 1.31e-03 4.42e+00  -5.7 6.65e-01  -3.7 1.00e+00 5.26e-01h  1
375  2.5548539e-02 2.62e-04 7.68e-01  -5.7 2.33e-01  -3.2 1.00e+00 9.70e-01f  1
376  3.1550045e-02 1.21e-03 6.07e+00  -5.7 1.30e+00  -3.7 6.77e-01 2.22e-01f  1
377  4.0770508e-02 9.00e-03 1.39e+01  -5.7 5.18e+02  -4.2 4.01e-04 5.02e-04f  1
378  5.6349896e-02 1.98e-02 1.17e+02  -5.7 4.82e-01  -3.8 3.23e-01 1.00e+00f  1
379  6.0645344e-02 1.95e-02 1.15e+02  -5.7 1.98e+01    -  7.16e-03 1.54e-02h  1
iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
380  6.6253489e-02 1.88e-02 1.11e+02  -5.7 5.16e+00  -4.2 3.69e-02 3.77e-02h  1
381  7.7560176e-02 1.77e-02 1.05e+02  -5.7 4.42e+00  -3.8 2.48e-02 5.62e-02f  1
382  9.8290149e-02 9.22e-03 5.73e+01  -5.7 5.59e-01  -3.4 3.41e-01 4.97e-01f  1
383  1.0750273e-01 8.91e-03 5.54e+01  -5.7 5.22e+00  -3.9 2.62e-02 3.38e-02h  1
384  1.2607729e-01 6.69e-03 4.26e+01  -5.7 1.21e+00  -3.4 1.24e-01 2.51e-01f  1
385  1.5119624e-01 5.88e-03 3.72e+01  -5.7 3.27e+00  -3.9 2.70e-02 1.53e-01f  1
386  1.5186515e-01 5.35e-03 3.50e+01  -5.7 1.85e+00    -  9.41e-02 1.00e-01h  1
387  1.5116854e-01 5.17e-03 3.34e+01  -5.7 3.65e+00    -  8.44e-02 3.31e-02h  1
388  1.5026803e-01 4.88e-03 3.08e+01  -5.7 2.86e+00    -  1.65e-01 5.59e-02h  1
389  1.4653227e-01 3.02e-03 1.66e+01  -5.7 7.02e-01    -  1.57e-01 4.11e-01h  1
iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
390  1.4732849e-01 2.26e-03 1.11e+01  -5.7 5.33e-01    -  3.58e-01 3.16e-01h  1
391  1.5068781e-01 7.33e-04 1.06e+01  -5.7 1.91e-01    -  8.26e-01 1.00e+00h  1
392  1.5212574e-01 3.18e-04 1.73e+00  -5.7 2.71e-01    -  8.46e-01 7.65e-01h  1
393  1.5206265e-01 1.44e-04 5.14e-01  -5.7 1.98e-01    -  1.00e+00 1.00e+00f  1
394  1.5115643e-01 3.49e-05 9.01e-02  -5.7 8.57e-02    -  1.00e+00 1.00e+00h  1
395  1.5088075e-01 1.95e-06 4.58e-03  -5.7 1.05e-02    -  1.00e+00 1.00e+00h  1
396  1.5086729e-01 5.32e-09 1.26e-05  -5.7 6.87e-04    -  1.00e+00 1.00e+00h  1
397  1.5102879e-01 5.96e-06 2.68e+01  -8.6 6.23e-02    -  8.51e-01 1.00e+00h  1
398  1.5105075e-01 1.22e-05 8.31e+00  -8.6 1.26e-01    -  6.90e-01 1.00e+00h  1
399  1.5105609e-01 8.66e-06 2.52e+00  -8.6 1.38e-01    -  6.97e-01 9.99e-01h  1
iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
400  1.5105708e-01 5.07e-06 9.56e-01  -8.6 2.65e-01    -  6.20e-01 4.93e-01h  1
401  1.5105818e-01 2.65e-06 5.51e-01  -8.6 3.28e-01    -  4.24e-01 7.56e-01f  1
402  1.5105850e-01 9.56e-07 2.77e-01  -8.6 6.52e-02    -  4.97e-01 7.78e-01f  1
403  1.5105872e-01 5.92e-07 2.12e-03  -8.6 9.67e-02    -  1.00e+00 1.00e+00f  1
404  1.5105875e-01 2.47e-07 1.58e-04  -8.6 3.97e-02    -  1.00e+00 1.00e+00h  1
405  1.5105875e-01 1.38e-09 4.24e-07  -8.6 4.29e-03    -  1.00e+00 1.00e+00h  1
406  1.5105875e-01 4.89e-13 7.06e-11  -8.6 8.47e-05    -  1.00e+00 1.00e+00h  1

Number of Iterations....: 406

(scaled)                 (unscaled)
Objective...............:  -1.5105874963825716e-01    1.5105874963825716e-01
Dual infeasibility......:   7.0607512892006952e-11    7.0607512892006952e-11
Constraint violation....:   4.8877568659122517e-13    4.8877568659122517e-13
Variable bound violation:   0.0000000000000000e+00    0.0000000000000000e+00
Complementarity.........:   2.5059039019713751e-09    2.5059039019713751e-09
Overall NLP error.......:   2.5059039019713751e-09    2.5059039019713751e-09

Number of objective function evaluations             = 1045
Number of objective gradient evaluations             = 306
Number of equality constraint evaluations            = 1045
Number of inequality constraint evaluations          = 1045
Number of equality constraint Jacobian evaluations   = 424
Number of inequality constraint Jacobian evaluations = 424
Number of Lagrangian Hessian evaluations             = 406
Total seconds in IPOPT                               = 27.202

EXIT: Optimal Solution Found.


Part of it is the numerics. Part of it is also how the du expressions are nested. Extracting the du into decision variables cuts the number of iterations required in half:

function integration_JuMP!(model, u, du, jc, ::TrapezoidalIntegrationRule)
@constraint(
model,
[t in 2:jc.n, k in 1:length(u)],
u[k][t] == u[k][t-1] + 0.5 * jc.Δt * (du[k][t-1] + du[k][t]),
)
return
end

function slow_sim(u0, jc::JumpConfig; Kd = 1e-6, k_on = 4e6, C_max = 1e-6)
n = jc.n
model = Model(Ipopt.Optimizer)
set_silent(model)
u = @variables(model,begin
0 <= φ_R[1:n] <= 0.45
0 <= c_pc[1:n]
0 <= C_bnd[1:n]
0 <= C_in[1:n]
0 <= M[1:n]
end);
for i in 1:length(u)
fix(u[i][1], u0[i]; force = true)
end
@variable(model, 0 <= α_R[1:n] <= 0.45)
@expressions(model, begin
C_Rf, 6.4e-5 .* φ_R .- C_bnd
φ_Rf, φ_R .- C_bnd ./ 6.4e-5
φ_M, 0.45 .- φ_R
γ, 10 .* c_pc ./ (c_pc .+ 0.03)
λ, γ .* φ_Rf
end)
@constraint(model, φ_Rf .>= 0)
k_off = k_on * Kd
du = [@variable(model, [1:n]) for _ in 1:length(u0)]
@constraints(model, begin
du[1] .== (α_R .- φ_R) .* λ
du[2] .== 5 .* (0.45 .- φ_R) .- (1 .+ c_pc) .* λ
du[3] .== k_on .* C_in .* C_Rf - k_off .* C_bnd .- λ .* C_bnd
du[4] .== 1.5 .* (C_max .- C_in) + k_off .* C_bnd .- k_on .* C_Rf .* C_in .- λ .* C_in
du[5] .== M .* λ
end)
integration_JuMP!(model, u, du, jc, TrapezoidalIntegrationRule())
@objective(model, Max, M[n])
optimize!(model)
@assert is_solved_and_feasible(model)
Rs = value.(φ_R)
a_Rs = value.(α_R)
cpcs = value.(c_pc)
C_bnds = value.(C_bnd)
C_ins = value.(C_in)
Ls = value.(λ)
return Rs, a_Rs, cpcs, C_bnds, C_ins, Ls
end

1 Like

thanks , that’s a huge improvement!

1 Like

what do you mean “decision variables”

what I got from seeing your changes is that lumping the expressions and the constraints into one macro call + not double defining du[i] with a second name reduces the overhead?

I created du as a decision variable (a column in the model).

This results in more variables in the model, but simpler constraints, since your integration constraints are now linear, and the nonlinearity is restricted to the @constraints block with du.

1 Like

You should carefully check whether the solutions make sense though. Ipopt has a tolerance of 1e-8, so you may be getting weird solutions if you have very large and very small data coefficients.

Perhaps a stupid question but why use JuMP to solve differential equations instead of DifferentialEquations.jl?

I am using JuMP to find the optimal control of the system – similarly to the rocket control eample in the docs.

In other scripts I am also using DifferentialEquations.jl to experiment with different control laws, but JuMP is giving me a reference for what is optimal

yesterday it ran much faster (~20s) with the changes you suggested, but today it’s back to taking much longer (~180s) and I haven’t changed anything.

here’s the Ipopt log:

This is Ipopt version 3.14.14, running with linear solver MUMPS 5.6.2.

Number of nonzeros in equality constraint Jacobian...:   109957
Number of nonzeros in inequality constraint Jacobian.:     4998
Number of nonzeros in Lagrangian Hessian.............:   112456

Total number of variables............................:    27495
variables with only lower bounds:     9996
variables with lower and upper bounds:     4999
variables with only upper bounds:        0
Total number of equality constraints.................:    24995
Total number of inequality constraints...............:     2500
inequality constraints with only lower bounds:     2500
inequality constraints with lower and upper bounds:        0
inequality constraints with only upper bounds:        0


I truncated the tail as before because it was simply far too long to fit

You must have changed something

The larger problem is that Ipopt assumes the problem is convex. It can struggle to solve nonconvex problems ike this one. You should expect that relatively minor changes to the formulation can have big changes to the solution speed.

You should also take a look at GitHub - infiniteopt/InfiniteOpt.jl: An intuitive modeling interface for infinite-dimensional optimization problems. which is intended for optimal control problems.

1 Like

ok, thanks!