 # Successive optimization in JuMP

Sorry for the long post but I thought it is better to give more information. The question itself is not too long! I have given the mathematical formulation below and also my Julia code.

Possible solution approach
If we ignore constraint (2) and suppose `adj` is a constant, then constraint 1 can be made linear using the big M method. Using candidate values for x, we can use a root finding algorithm to find adj. We then update the value of `adj` in constraint (1), optimize again and use candidate values of x to find `adj`. The process is repeated until changes in `adj` are less than a certain threshold or we hit the limit for maximum number of iterations.
For this to work, I provide an initial guess value for adj. I have noticed that the final solution is somewhat sensitive to this guess value. I would like advice on the following two parts:

• Part1: In the function, `successive_linear` I define an array ` adj_trials = []` (as I don’t know its dimension) to hold values for `adj` and I use `push!` to store values. A while loop is used to iterate with the condition for changes in `adj` to be within certain tolerance and number of iterations under a limit. Is there a better way of using the while loop with these conditions, either the way I store values in the array `adj_trails` or how I test for the termination condition? Should I also add changes in value of the objective function as an extra termination condition, i.e. if changes are less than a threshold?

• Part2: Any suggestions for a better formulation of this problem or how the final solution could be made less sensitive to the initial? I don’t know if it is the Cbc solver that gives non-deterministic results or some other factor is playing up.

``````using JuMP, Cbc, Roots

function some_data()
cost = [500,400,300,600,450]
A= [100 120 130 140 150; 100 150 200 300 90; 100 200 150 140 150; 100 120 130 140 150; 200 300 400 500 200]
b = [200, 300, 400, 1000, 500, 20, 30, 3]
if size(A,2) < length(b)
A = hcat(A,zeros(size(A,1),length(b)-size(A,2)))
elseif length(b) < size(A,2)
b = vcat(b,zeros(size(A,2)-length(b)))
end
Factor1 = [0.001, 0.015, 0.02, 0.025, 0.028, 0.03,0.031, 0.035]
Factor2 = [0.0015, 0.017, 0.022, 0.024,0.03, 0.0301, 0.0306,0.033]
t = length(b)
n = size(A,1)

return cost, A, b, Factor1, Factor2, n , t
end

model = Model(Cbc.Optimizer)

@variables(model, begin
0 <= x[1:n] <= 1
y[1:t]
y_p[1:t] >= 0
y_n[1:t] >= 0
z[1:t], Bin
end)

@expression(model, expr[i=1:t], A[:, i]' * x - b[i] + y[i])
fix(y, 0.0)

for i in 2:t
@constraints(model, begin
y_p[i] - y_n[i] == expr[i - 1]
y[i] == (1+Factor1[i-1]) * y_p[i] - (1+Factor2[i-1]+adj_factor_trial) * y_n[i]
y_p[i] <= M * z[i]
y_n[i] <= M * (1 - z[i])
end)
end
@constraint(model, con_zero, expr[t] == 0)
@objective(model,Min,sum(cost[i]*x[i] for i in 1:n))
return x,y,y_p,y_n,z, model
end

f(x) = sum(c[i]/(1+Factor2[i]+x)^i for i in 1:length(c))-cost
end

cost, A, b, Factor1, Factor2,n , t = some_data()
M=1000000
x,y,y_p,y_n,z, model = indicator_cons(cost,A,b,Factor1,Factor2,M, n, t,adj_guess)
optimize!(model)
x_r = value.(x)
y_r = value.(y)
y_p_r= value.(y_p)
y_n_r = value.(y_n)
z_r = value.(z)

iter = 1
println("Iteration \$(iter)")
while true
iter += 1
x,y,y_p,y_n,z, model = indicator_cons(cost,A,b,Factor1,Factor2,M, n, t,adj)
JuMP.set_start_value.(x, x_r)
JuMP.set_start_value.(y, y_r)
JuMP.set_start_value.(y_p, y_p_r)
JuMP.set_start_value.(y_n, y_n_r)
JuMP.set_start_value.(z, z_r)

optimize!(model)
x_r = value.(x)
y_r = value.(y)
y_p_r= value.(y_p)
y_n_r = value.(y_n)
z_r = value.(z)

println("Iteration \$(iter)")
println("Tolerance \$(tol)")

if abs(tol) <= ftol || iter > maxiter
break
end
end

end

``````

I thought I recognized this problem:

Is there a better way of using the while loop with these conditions

A while loop like this is fine. You could also define

``````iter = 1
tol = Inf
while abs(tol) > ftol && iter <= maxiter
iter += 1
# ...
end
``````

how the final solution could be made less sensitive to the initial?

I haven’t run your code, but my guess is that there are lots of local optima. Are you finding different `x` solutions with the same objective values? Or are the objective values different as well?

What does `adj` represent in reality?

Thank you.

Sorry, I used the same example as in my previous post for a different question, but didn’t referenced it.

I find different x values with different objective values, so objective values are different as well.
Actually `adj` is what I want to maximize, but I thought it will be harder to do that as (2) is nonconvex so I came up with this objective function which is linear. `Adj` represents a flat addition to some parameters `Factor2`, but it depends upon x values. A higher value of `adj` is better.

If I define the `adj` as JuMP variable `@variable(model,adj)` and change the following line of code ` y[i] == (1+Factor1[i-1]) * y_p[i] - (1+Factor2[i-1]+adj) * y_n[i]`. This will make this constraint quadratic so it will no be supported by Cbc. Perhaps product of `adj` and `y_n` can be linearized. Then changing the objective function to maximize `adj`, we can find upper bound on value of `adj`. Not sure if it will help with the solving this problem as I am ignoiring constraint (2).

Can we write a JuMP constraint with if then else control statements and let a nonlinear solver to do the magic? I tried it but JuMP doesnt seem to support it, but I could be getting the syntax wrong.

Gurobi lets you have bilinear terms like this.

But the `^(1/12)` is a problem you won’t be able to easily work-around.

I find different x values with different objective values, so objective values are different as well.

Did you try plotting the objective value against `adj`? What does that look like?

If `adj` is a scalar, and bigger values are better, you could just do a line search.

Thank you

Do you mean giving different guess values for `adj` and then running the optimization, and then plotting the ‘optimal’ `adj` and objective value? I did a plot of adj and objective based on that. It shows a closed arc type of shape at first sight, but if I examine it carefully the difference in results based on different guess values of `adj` is relatively small, but there are differences in x values. I don’t know if this is due to multiple local optima, but more likely looks like sensitivity of the solver. The data I am using has small and large values, which could be causing numerical issues.

Yes. If it is an arc, you could use something like golden-section search: Line search - Wikipedia

The data I am using has small and large values, which could be causing numerical issues.

Try a different solver like Gurobi and check you get similar results. Gurobi is more likely to warn you if it encounters numerical issues.

Thanks. I will explore these line search methods and see if it works!