# Successive optimization in JuMP
**URL:** https://discourse.julialang.org/t/successive-optimization-in-jump/64056
**Category:** Optimization (Mathematical)
**Tags:** jump, optimization
**Created:** [July 4, 2021, 11:41pm UTC](https://discourse.julialang.org/t/successive-optimization-in-jump/64056 "2021-07-04T23:41:24Z")
**Posts on this page:** 7
**Page:** 1
### Author:  [@Maria95](https://discourse.julialang.org/u/Maria95)
#### Post date: [July 4, 2021, 11:41pm UTC](https://discourse.julialang.org/t/successive-optimization-in-jump/64056/1 "2021-07-04T23:41:24Z")
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.
```julia
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
function indicator_cons(cost,A,b,Factor1,Factor2,M, n, t,adj_guess)
model = Model(Cbc.Optimizer)
adj_factor_trial = -1+(1+adj_guess)^(1/12)
@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[1], 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
function factor_adj(cost,c,Factor2)
f(x) = sum(c[i]/(1+Factor2[i]+x)^i for i in 1:length(c))-cost
adj = find_zero(f, (0,1), Bisection())
return adj
end
function successive_linear(;ftol=1.0e-9,adj_guess = 0.4,maxiter=20)
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)
adj = factor_adj(cost'*x_r,A'*x_r,Factor2)
adj_trials = []
push!(adj_trials,adj)
iter = 1
println("Iteration $(iter)")
println("Adjustment factor $(adj)")
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)
adj = factor_adj(cost'*x_r,A'*x_r,Factor2)
push!(adj_trials,adj)
println("Iteration $(iter)")
println("Adjustment factor $(adj)")
tol = adj_trials[iter] - adj_trials[iter-1]
println("Tolerance $(tol)")
if abs(tol) <= ftol || iter > maxiter
break
end
end
return x_r, adj_trials
end
x_r, adj_trials = successive_linear(;ftol=1.0e-9,adj_guess = 0.01,maxiter=20)
```
---
### Author:  [@odow](https://discourse.julialang.org/u/odow)
#### Post date: [July 5, 2021, 12:20am UTC](https://discourse.julialang.org/t/successive-optimization-in-jump/64056/2 "2021-07-05T00:20:43Z")
I thought I recognized this problem:
> [@Updating expressions for indicator constraints in JuMP](https://discourse.julialang.org/t/updating-expressions-for-indicator-constraints-in-jump/59059):
>
> Hello, I am trying to implement a constraint with a conditional statement that I am making piecewise affine by introducing binary variables and using the big M method. I know that JuMP supports indicator constraints but I have not been able to find any relevant examples for more complicated use cases such as this one. My question may seem long, but I wanted to give enough information so that it is clear. When I run my code and look at how the expressions are being built, it does not look like t…
> Is there a better way of using the while loop with these conditions
A while loop like this is fine. You could also define
```Julia
iter = 1
tol = Inf
while abs(tol) > ftol && iter <= maxiter
iter += 1
# ...
tol = adj_trials[iter] - adj_trials[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?
---
### Author:  [@Maria95](https://discourse.julialang.org/u/Maria95)
#### Post date: [July 5, 2021, 10:36pm UTC](https://discourse.julialang.org/t/successive-optimization-in-jump/64056/3 "2021-07-05T22:36:42Z")
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.
---
### Author:  [@odow](https://discourse.julialang.org/u/odow)
#### Post date: [July 5, 2021, 11:54pm UTC](https://discourse.julialang.org/t/successive-optimization-in-jump/64056/4 "2021-07-05T23:54:24Z")
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.
---
### Author:  [@Maria95](https://discourse.julialang.org/u/Maria95)
#### Post date: [July 6, 2021, 10:43pm UTC](https://discourse.julialang.org/t/successive-optimization-in-jump/64056/5 "2021-07-06T22:43:26Z")
Thank you
> [@odow](#):
>
> Did you try plotting the objective value against `adj` ?
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.
---
### Author:  [@odow](https://discourse.julialang.org/u/odow)
#### Post date: [July 6, 2021, 11:01pm UTC](https://discourse.julialang.org/t/successive-optimization-in-jump/64056/6 "2021-07-06T23:01:13Z")
> [@Maria95](#):
>
> Do you mean giving different guess values for `adj` and then running the optimization, and then plotting the ‘optimal’ `adj` and objective value?
Yes. If it is an arc, you could use something like golden-section search: [Line search - Wikipedia](https://en.wikipedia.org/wiki/Line_search)
> 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.
---
### Author:  [@Maria95](https://discourse.julialang.org/u/Maria95)
#### Post date: [July 7, 2021, 11:03pm UTC](https://discourse.julialang.org/t/successive-optimization-in-jump/64056/7 "2021-07-07T23:03:43Z")
Thanks. I will explore these line search methods and see if it works!