# Matrix inversion inside JuMP objective

Hello,

I am trying to solve the following problem using JuMP: find a set of values `u` subject to different linear constraints (namely positivity and mass conservation) such that the distance to a know solution of a linear system is minimal: \min ||x_{eq} + A(u)^{-1}b(u)||_2
However, I’ve tried several things but I am still unable to implement the inverse in the objective function without error.

Here is a chunk of my code:

``````using JuMP, LinearAlgebra
using Ipopt: Ipopt

age_model = Model(Ipopt.Optimizer)

vol = rand(9)
xeq = rand(5)

xinf = rand(4)

function exchange(u)
[
0 u u u 0 0 0 0 0;
u 0 0 0 u 0 0 0 0;
u 0 0 u 0 u 0 u 0;
u 0 u 0 u 0 u u 0;
0 u 0 u 0 0 0 0 u;
0 0 u 0 0 0 u u 0;
0 0 0 u 0 u 0 u 0;
0 0 u u 0 u u 0 u;
0 0 0 0 u 0 0 u 0
]
end

function A(u)
[
-((u + u) / vol) 0 u/vol 0 0;
0 -((u + u + u + u + u) / vol) u/vol u/vol 0;
u/vol u/vol -((u + u + u) / vol) 0 u/vol;
0 u/vol 0 -((u + u + u) / vol) 0;
0 0 u/vol 0 -((u + u) / vol)
]
end

function b(u)
[
(u * xinf) / vol,
(u * xinf + u * xinf + u * xinf) / vol,
0,
(u * xinf + u * xinf) / vol,
(u * xinf) / vol,
]
end

function objfun(u)
norm(xeq .+ A(u)\b(u))
end

@variable(age_model, u[1:30] >= 0)

@constraint(age_model, vec(sum(exchange(u), dims = 1)) .== vec(sum(exchange(u), dims = 2)))

register(age_model, :objfun, 1, objfun; autodiff=true)

@NLobjective(age_model, Min, objfun(u...))

optimize!(age_model)
``````

Thanks,
L.

The good news is that any term of the type v = A^{-1} b can be written as a constraint A v = b with a new variable v. So roughly:

``````@variable(model, v[1:n])

@NLconstraint(model, A v .== b)
``````

for a matrix `A` with `n` columns. You’ll likely need to expand the last line to work with the scalar nonlinear expressions that involve variables in the coefficients.

4 Likes

Adding a variable did the trick, thanks !

Another (maybe not really related question), do you know if there is an ‘optimized’ method to launch many optimizations, each time with a perturbed starting point?

I.e: I have `start = u0[i])` when defining the variables, but is there a way to set it to `u0[i] + 0.1*rand()` without having to redefine the model in a loop?

You could probably use a Parameter:

You don’t need to use a parameter, just set a new start value:

``````for i in 1:30
set_start_value(u[i], u0[i] + 0.1 * rand())
end
optimize!(age_model)
``````

Docs: Variables · JuMP

2 Likes