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_linearI define an array
adj_trials = (as I don’t know its dimension) to hold values for
adjand I use
push!to store values. A while loop is used to iterate with the condition for changes in
adjto 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_trailsor 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 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, 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)