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 arrayadj_trials = [](as I don’t know its dimension) to hold values foradjand I usepush!to store values. A while loop is used to iterate with the condition for changes inadjto 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 arrayadj_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[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)