Matrix inversion inside JuMP objective

Specific Domains Optimization (Mathematical)
1

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[1] u[2] u[3] 0 0 0 0 0;
		u[4] 0 0 0 u[5] 0 0 0 0;
		u[6] 0 0 u[7] 0 u[8] 0 u[9] 0;
		u[10] 0 u[11] 0 u[12] 0 u[13] u[14] 0;
		0 u[15] 0 u[16] 0 0 0 0 u[17];
		0 0 u[18] 0 0 0 u[19] u[20] 0;
		0 0 0 u[21] 0 u[22] 0 u[23] 0;
		0 0 u[24] u[25] 0 u[26] u[27] 0 u[28];
		0 0 0 0 u[29] 0 0 u[30] 0
	]
end

function A(u)
	[
		-((u[1] + u[15]) / vol[2]) 0 u[5]/vol[2] 0 0;
		0 -((u[3] + u[7] + u[16] + u[21] + u[25]) / vol[4]) u[12]/vol[4] u[13]/vol[4] 0;
		u[15]/vol[5] u[16]/vol[5] -((u[5] + u[12] + u[29]) / vol[5]) 0 u[17]/vol[5];
		0 u[21]/vol[7] 0 -((u[13] + u[19] + u[27]) / vol[7]) 0;
		0 0 u[29]/vol[9] 0 -((u[17] + u[28]) / vol[9])
	]
end

function b(u)
	[
		(u[4] * xinf[1]) / vol[2],
		(u[10] * xinf[1] + u[11] * xinf[2] + u[14] * xinf[4]) / vol[4],
		0,
		(u[22] * xinf[3] + u[23] * xinf[4]) / vol[7],
		(u[30] * xinf[4]) / vol[9],
	]
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.

2

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
4

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?

5

You could probably use a Parameter:

6

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