I am trying to solve a problem where I wish to parameterize certain constants in my model for use in DiffOpt.jl. I use the constants to recursively generate a signal that I use elsewhere in the model. Formulating this problem using JuMP parameters takes much, much longer than if I were to just calculate the signal outside of an expression and formulate the optimization problem that way.
Here is a minimal example.
ϕ0=[1.6738595158677403, -0.4498714736488657, -0.225851131510401])
init = [0.1, 0.2, 0.3];
horz = 50;
# create model
model = Model(optimizer_with_attributes(Ipopt.Optimizer))
# create variables
@variables(model, begin
x[t=1:horz]
ϕ[i=1:3] in Parameter(ϕ0[i])
end)
A = [ϕ[1] ϕ[2] ϕ[3];
1 0 0;
0 1 0];
@expression(model, Pref[n=1:3, t=1:horz],
LinearAlgebra.dot( (A^(t-1))[n,:], init)
);
@time @objective(model, Min, sum((Pref[1, t] - x[t]).^2 for t=1:horz))
@time optimize!(model)
assert_is_solved_and_feasible(model)
The bulk of the runtime is spent on the objective macro and the optimize, although the vast majority of the time in the optimize run is spent outside of Ipopt. This does speed up on subsequent runs of the model, but I would like it to be faster, if possible.
6.254010 seconds (59.23 M allocations: 1.961 GiB, 32.15% gc time, 1.09% compilation time)
[...]
Total seconds in IPOPT = 0.206
EXIT: Optimal Solution Found.
6.943605 seconds (50.29 M allocations: 2.534 GiB, 42.47% gc time)
Are there any places where I can speed up the model construction and solving?
By the time you get to A^49, the expressions are veeeeery large.
You can work around this by adding intermediate variables instead of a large nested expression. Assuming I’ve understood your example correctly, something like this:
using JuMP, Ipopt
ϕ0 = [1.6738595158677403, -0.4498714736488657, -0.225851131510401]
init = [0.1, 0.2, 0.3]
horz = 50
model = Model(Ipopt.Optimizer)
@variables(model, begin
x[t in 1:horz]
ϕ[i in 1:3] in Parameter(ϕ0[i])
pref[1:3, 1:horz]
end)
A = [ϕ[1] ϕ[2] ϕ[3]; 1 0 0; 0 1 0];
@constraints(model, begin
pref[:, 1] .== init
[t in 2:horz], pref[:, t] .== A * pref[:, t - 1]
end)
@objective(model, Min, sum((pref[1, t] - x[t])^2 for t in 1:horz))
optimize!(model)
assert_is_solved_and_feasible(model)