Skip to content

Commit

Permalink
update
Browse files Browse the repository at this point in the history
  • Loading branch information
@jiedxu
jiedxu committed Aug 12, 2019
1 parent 87c6b78 commit e17952e
Show file tree
Hide file tree
Showing 8 changed files with 89 additions and 89 deletions.
36 changes: 18 additions & 18 deletions src/benders/lshaped.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -20,7 +20,7 @@ function lshaped(; n_x, n_y, vec_min_y, vec_max_y, vec_f,
# Define Master problem
n_constraint = length(mat3_w[1, :, 1])
num_s = length(mat3_t[:, 1, 1])
model_mas = Model(solver = GLPKSolverMIP())
model_mas = Model(with_optimizer(GLPK.Optimizer))
@variable(model_mas, q)
@variable(model_mas, vec_y[1: n_y], Int)
@objective(model_mas, Min, (transpose(vec_f) * vec_y + q)[1])
Expand All @@ -35,9 +35,9 @@ function lshaped(; n_x, n_y, vec_min_y, vec_max_y, vec_f,
@constraint(model_mas, (mat_e1 * vec_y)[1] >= e2)
end
@constraint(model_mas, (mat_e1 * vec_y)[1] >= - q + e2)
solve(model_mas)
vec_result_y = getvalue(vec_y)
return getobjectivevalue(model_mas)
optimize!(model_mas)
vec_result_y = value(vec_y)
return objective_value(model_mas)
end


Expand All @@ -46,17 +46,17 @@ function lshaped(; n_x, n_y, vec_min_y, vec_max_y, vec_f,
@variable(model_sub, vec_u[1: n_constraint] >= 0)
@objective(model_sub, Max, (transpose(vec_h - mat_t * vec_yBar) * vec_u)[1])
constraintsForDual = @constraint(model_sub, transpose(mat_w) * vec_u .<= vec_c)
solution_sub = solve(model_sub)
solution_sub = optimize!(model_sub)
print("------------------------------ Sub Problem ------------------------------\n") # , model_sub)
vec_uBar = getvalue(vec_u)
vec_uBar = value(vec_u)
if solution_sub == :Optimal
vec_result_x = zeros(length(vec_c))
vec_result_x = getdual(constraintsForDual)
return (true, getobjectivevalue(model_sub), vec_uBar, vec_result_x)
vec_result_x = dual(constraintsForDual)
return (true, objective_value(model_sub), vec_uBar, vec_result_x)
end
if solution_sub == :Unbounded
print("Not solved to optimality because feasible set is unbounded.\n")
return (false, getobjectivevalue(model_sub), vec_uBar, repeat([NaN], length(vec_c)))
return (false, objective_value(model_sub), vec_uBar, repeat([NaN], length(vec_c)))
end
if solution_sub == :Infeasible
print("Not solved to optimality because infeasibility. Something is wrong.\n")
Expand All @@ -71,10 +71,10 @@ function lshaped(; n_x, n_y, vec_min_y, vec_max_y, vec_f,
@objective(model_ray, Max, 1)
@constraint(model_ray, (transpose(vec_h - mat_t * vec_yBar) * vec_u)[1] == 1)
@constraint(model_ray, transpose(mat_w) * vec_u .<= 0)
solve(model_ray)
optimize!(model_ray)
print("------------------------------ Ray Problem ------------------------------\n") # , model_ray)
vec_uBar = getvalue(vec_u)
obj_ray = getobjectivevalue(model_ray)
vec_uBar = value(vec_u)
obj_ray = objective_value(model_ray)
return (obj_ray, vec_uBar)
end

Expand Down Expand Up @@ -119,15 +119,15 @@ function lshaped(; n_x, n_y, vec_min_y, vec_max_y, vec_f,
mat_e1 = sum(vec_pi[s] * (transpose(mat_uBar[s, :]) * mat3_t[s, :, :])[1] for s = 1: num_s)
e2 = sum(vec_pi[s] * (transpose(mat_uBar[s, :]) * mat_h[s, :])[1] for s = 1: num_s)
obj_mas = solve_master(mat_e1, e2, vec_bool_solutionSubModel[1])
vec_yBar = getvalue(vec_y)
vec_yBar = value(vec_y)
# 3. Compare the bounds and decide whether to stop
boundLow = max(boundLow, obj_mas)
dict_boundUp[timesIteration] = boundUp
dict_boundLow[timesIteration] = boundLow
if vec_bool_solutionSubModel[1]
dict_obj_mas[timesIteration] = obj_mas
dict_obj_sub[timesIteration] = obj_sub
result_q = getvalue(q)
result_q = value(q)
dict_q[timesIteration] = result_q
println("------------------ Result in $(timesIteration)-th Iteration with Sub ",
"-------------------\n", "boundUp: $(round(boundUp, digits = 5)), ",
Expand All @@ -136,7 +136,7 @@ function lshaped(; n_x, n_y, vec_min_y, vec_max_y, vec_f,
else
dict_obj_mas[timesIteration] = obj_mas
dict_obj_ray[timesIteration] = obj_sub
result_q = getvalue(q)
result_q = value(q)
dict_q[timesIteration] = result_q
println("------------------ Result in $(timesIteration)-th Iteration with Ray ",
"-------------------\n", "boundUp: $(round(boundUp, digits = 5)), ",
Expand All @@ -145,7 +145,7 @@ function lshaped(; n_x, n_y, vec_min_y, vec_max_y, vec_f,
end
timesIteration += 1
end # -----------------------------------------------------------------------------------------------------
println("obj_mas: $(getobjectivevalue(model_mas))")
println("obj_mas: $(objective_value(model_mas))")
println("----------------------------- Master Problem ----------------------------\n")
# println(model_mas)
println("-------------------------------------------------------------------------\n",
Expand All @@ -154,8 +154,8 @@ function lshaped(; n_x, n_y, vec_min_y, vec_max_y, vec_f,
println("boundUp: $(round(boundUp, digits = 5)), boundLow: $(round(boundLow, digits = 5)), ",
"difference: $(round(boundUp - boundLow, digits = 5))")
println("vec_x: $vec_result_x")
vec_result_y = getvalue(vec_y)
result_q = getvalue(q)
vec_result_y = value(vec_y)
result_q = value(q)
println("vec_y: $vec_result_y")
println("result_q: $result_q")
println("-------------------------------------------------------------------------\n",
Expand Down
36 changes: 18 additions & 18 deletions src/benders/milp.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -13,7 +13,7 @@ function milp(; n_x, n_y, vec_min_y, vec_max_y, vec_c, vec_f, vec_b, mat_a, mat_
"-------------------------------------------------------------------------\n")
# Define Master problem
n_constraint = length(mat_a[:, 1])
model_mas = Model(solver = GLPKSolverMIP())
model_mas = Model(with_optimizer(GLPK.Optimizer))
@variable(model_mas, q)
@variable(model_mas, vec_y[1: n_y], Int)
@objective(model_mas, Min, (transpose(vec_f) * vec_y + q)[1])
Expand All @@ -28,9 +28,9 @@ function milp(; n_x, n_y, vec_min_y, vec_max_y, vec_c, vec_f, vec_b, mat_a, mat_
@constraint(model_mas, (transpose(vec_uBar) * (vec_b - mat_b * vec_y))[1] <= 0)
end
@constraint(model_mas, (transpose(vec_uBar) * (vec_b - mat_b * vec_y))[1] <= q)
solve(model_mas)
vec_result_y = getvalue(vec_y)
return getobjectivevalue(model_mas)
optimize!(model_mas)
vec_result_y = value(vec_y)
return objective_value(model_mas)
end


Expand All @@ -39,17 +39,17 @@ function milp(; n_x, n_y, vec_min_y, vec_max_y, vec_c, vec_f, vec_b, mat_a, mat_
@variable(model_sub, vec_u[1: n_constraint] >= 0)
@objective(model_sub, Max, (transpose(vec_b - mat_b * vec_yBar) * vec_u)[1])
constraintsForDual = @constraint(model_sub, transpose(mat_a) * vec_u .<= vec_c)
solution_sub = solve(model_sub)
solution_sub = optimize!(model_sub)
print("------------------------------ Sub Problem ------------------------------\n") # , model_sub)
vec_uBar = getvalue(vec_u)
vec_uBar = value(vec_u)
if solution_sub == :Optimal
vec_result_x = zeros(length(vec_c))
vec_result_x = getdual(constraintsForDual)
return (true, getobjectivevalue(model_sub), vec_uBar, vec_result_x)
vec_result_x = dual(constraintsForDual)
return (true, objective_value(model_sub), vec_uBar, vec_result_x)
end
if solution_sub == :Unbounded
print("Not solved to optimality because feasible set is unbounded.\n")
return (false, getobjectivevalue(model_sub), vec_uBar, repeat([NaN], length(vec_c)))
return (false, objective_value(model_sub), vec_uBar, repeat([NaN], length(vec_c)))
end
if solution_sub == :Infeasible
print("Not solved to optimality because infeasibility. Something is wrong.\n")
Expand All @@ -64,10 +64,10 @@ function milp(; n_x, n_y, vec_min_y, vec_max_y, vec_c, vec_f, vec_b, mat_a, mat_
@objective(model_ray, Max, 1)
@constraint(model_ray, (transpose(vec_b - mat_b * vec_yBar) * vec_u)[1] == 1)
@constraint(model_ray, transpose(mat_a) * vec_u .<= 0)
solve(model_ray)
optimize!(model_ray)
print("------------------------------ Ray Problem ------------------------------\n") # , model_ray)
vec_uBar = getvalue(vec_u)
obj_ray = getobjectivevalue(model_ray)
vec_uBar = value(vec_u)
obj_ray = objective_value(model_ray)
return (obj_ray, vec_uBar)
end

Expand Down Expand Up @@ -101,14 +101,14 @@ function milp(; n_x, n_y, vec_min_y, vec_max_y, vec_c, vec_f, vec_b, mat_a, mat_
(obj_ray, vec_uBar) = solve_ray(vec_yBar, n_constraint, vec_b, mat_b, mat_a)
end
obj_mas = solve_master(vec_uBar, bool_solutionSubModel)
vec_yBar = getvalue(vec_y)
vec_yBar = value(vec_y)
boundLow = max(boundLow, obj_mas)
dict_boundUp[timesIteration] = boundUp
dict_boundLow[timesIteration] = boundLow
if bool_solutionSubModel
dict_obj_mas[timesIteration] = obj_mas
dict_obj_sub[timesIteration] = obj_sub
result_q = getvalue(q)
result_q = value(q)
dict_q[timesIteration] = result_q
println("------------------ Result in $(timesIteration)-th Iteration with Sub ",
"-------------------\n", "boundUp: $(round(boundUp, digits = 5)), ",
Expand All @@ -117,7 +117,7 @@ function milp(; n_x, n_y, vec_min_y, vec_max_y, vec_c, vec_f, vec_b, mat_a, mat_
else
dict_obj_mas[timesIteration] = obj_mas
dict_obj_ray[timesIteration] = obj_ray
result_q = getvalue(q)
result_q = value(q)
dict_q[timesIteration] = result_q
println("------------------ Result in $(timesIteration)-th Iteration with Ray ",
"-------------------\n", "boundUp: $(round(boundUp, digits = 5)), ",
Expand All @@ -126,7 +126,7 @@ function milp(; n_x, n_y, vec_min_y, vec_max_y, vec_c, vec_f, vec_b, mat_a, mat_
end
timesIteration += 1
end
println("obj_mas: $(getobjectivevalue(model_mas))")
println("obj_mas: $(objective_value(model_mas))")
println("----------------------------- Master Problem ----------------------------\n")
# println(model_mas)
println("-------------------------------------------------------------------------\n",
Expand All @@ -135,8 +135,8 @@ function milp(; n_x, n_y, vec_min_y, vec_max_y, vec_c, vec_f, vec_b, mat_a, mat_
println("boundUp: $(round(boundUp, digits = 5)), boundLow: $(round(boundLow, digits = 5)), ",
"difference: $(round(boundUp - boundLow, digits = 5))")
println("vec_x: $vec_result_x")
vec_result_y = getvalue(vec_y)
result_q = getvalue(q)
vec_result_y = value(vec_y)
result_q = value(q)
println("vec_y: $vec_result_y")
println("result_q: $result_q")
println("-------------------------------------------------------------------------\n",
Expand Down
10 changes: 5 additions & 5 deletions src/dw/mas.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -57,7 +57,7 @@ function addColToMas(modMas::JuMP.Model, modSub::ModelSub, vec_x_result, vecCons
vecConsTouched = ConstraintRef[]
vals = Float64[]
for i = 1: m
# only insert non-zero elements (this saves memory and may make the modMas problem easier to solve)
# only insert non-zero elements (this saves memory and may make the modMas problem easier to optimize!)
if A0x[i] != 0
push!(vecConsTouched, vecConsRef[i])
push!(vals, A0x[i])
Expand All @@ -80,13 +80,13 @@ end


function solveMas(modMas, vecConsRef, vecConsConvex)
status = solve(modMas)
status = optimize!(modMas)
if status != :Optimal
throw("Error: Non-optimal modMas-problem status")
end
vecPi = getdual(vecConsRef)
vecKappa = getdual(vecConsConvex)
obj_master = getobjectivevalue(modMas)
vecPi = dual(vecConsRef)
vecKappa = dual(vecConsConvex)
obj_master = objective_value(modMas)
## ensure that vecPi and vecKappa are row vectors
vecPi = reshape(vecPi, 1, length(vecPi))
vecKappa = reshape(vecKappa, 1, length(vecKappa))
Expand Down
4 changes: 2 additions & 2 deletions src/dw/optim.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -70,8 +70,8 @@ function doOptim(vecModelSub, modMas, vecConsRef, vecConsConvex, vecLambda, dual
end
## 5, Print the Result
println("Optimization Done After $(iter-1) Iterations.\n",
"objMaster = $(getobjectivevalue(modMas))")
vecLambdaResult = getvalue(vecLambda)
"objMaster = $(objective_value(modMas))")
vecLambdaResult = value(vecLambda)
return (vecLambdaResult, extremePointForSub, extremePoints)
end

Expand Down
6 changes: 3 additions & 3 deletions src/dw/sub.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -74,11 +74,11 @@ function solveSub(modSub::ModelSub, vecPi, kappa)
piA0 = vecPi * modSub.mat_e
@objective(modSub.mod, Max, sum(modSub.vec_l[j] * modSub.vec_x[j] for j = 1: n) -
sum(piA0[1, j] * modSub.vec_x[j] for j = 1: n) - kappa)
status = solve(modSub.mod)
status = optimize!(modSub.mod)
if status != :Optimal
throw("Error: Non-optimal sub-problem status")
end
costReduce = getobjectivevalue(modSub.mod)
vec_x_result = getvalue(modSub.vec_x)
costReduce = objective_value(modSub.mod)
vec_x_result = value(modSub.vec_x)
return (costReduce, vec_x_result)
end
16 changes: 8 additions & 8 deletions src/milp/func.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -95,10 +95,10 @@ function solveLinear(vec_c, vec_b, mat_aCap)
@variable(model, vec_x[1: n_x] >= 0)
@objective(model, Min, (transpose(vec_c) * vec_x)[1])
constraintsForDual = @constraint(model, mat_aCap * vec_x .>= vec_b)
solve(model)
vec_result_u = getdual(constraintsForDual)
obj = getobjectivevalue(model)
vec_result_x = getvalue(vec_x)
optimize!(model)
vec_result_u = dual(constraintsForDual)
obj = objective_value(model)
vec_result_x = value(vec_x)
return obj, vec_result_x, vec_result_u
end

Expand All @@ -108,9 +108,9 @@ function solveMix(n_x, vec_c, vec_b, mat_aCap, vec_f, mat_b)
@variable(model, vec_x[1: n_x] >= 0)
@objective(model, Min, (transpose(vec_c) * vec_x)[1])
constraintsForDual = @constraint(model, mat_aCap * vec_x .>= vec_b)
solve(model)
vec_result_u = getdual(constraintsForDual)
obj = getobjectivevalue(model)
vec_result_x = getvalue(vec_x)
optimize!(model)
vec_result_u = dual(constraintsForDual)
obj = objective_value(model)
vec_result_x = value(vec_x)
return obj, vec_result_x, vec_result_u
end
30 changes: 15 additions & 15 deletions src/robust/box.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -19,7 +19,7 @@ function milpBoxBudget(; num_x, num_y, vec_min_y, vec_max_y, vec_c, vec_f, vec_b
"---------------------------------------------------------\n")
println("Input: vec_gammaCap = $vec_gammaCap.")
#
model = Model(solver = GLPKSolverMIP())
model = Model(with_optimizer(GLPK.Optimizer))
@variable(model, vec_y[1: num_y], Int)
@variable(model, vec_x[1: num_x] >= 0)
@objective(model, Min, (transpose(vec_c) * vec_x + transpose(vec_f) * vec_y)[1])
Expand All @@ -38,15 +38,15 @@ function milpBoxBudget(; num_x, num_y, vec_min_y, vec_max_y, vec_c, vec_f, vec_b
sum(mat_mu[i, :]) + transpose(mat_b_bar[i, :]) * vec_y <= vec_b_bar[i])
@constraint(model, vec_lambda[i] .+ hcat(mat_mu[i, :]) .>= hcat(mat_a_hat[i, :]) .* vec_z[:])
end
solve(model)
obj_result = getobjectivevalue(model)
vec_result_y = getvalue(vec_y)
optimize!(model)
obj_result = objective_value(model)
vec_result_y = value(vec_y)
println("Result: vec_y = $vec_result_y.")
vec_result_x = getvalue(vec_x)
vec_result_x = value(vec_x)
println("Result: vec_x = $vec_result_x.")
vec_result_z = getvalue(vec_z)
mat_result_mu = getvalue(mat_mu)
vec_result_lambda = getvalue(vec_lambda)
vec_result_z = value(vec_z)
mat_result_mu = value(mat_mu)
vec_result_lambda = value(vec_lambda)
println("---------------------------------------------------------\n",
" 2/2. Nominal Ending\n",
"---------------------------------------------------------\n")
Expand All @@ -69,7 +69,7 @@ function lpBoxBudget(; num_x, vec_c, vec_b, mat_a, mat_a_bar, mat_a_hat, vec_b_b
"---------------------------------------------------------\n")
println("Input: vec_gammaCap = $vec_gammaCap.")
#
model = Model(solver = GLPKSolverMIP())
model = Model(with_optimizer(GLPK.Optimizer))
@variable(model, vec_x[1: num_x] >= 0)
@objective(model, Min, (transpose(vec_c) * vec_x)[1])
# Transformation of Box Uncertainty.
Expand All @@ -84,13 +84,13 @@ function lpBoxBudget(; num_x, vec_c, vec_b, mat_a, mat_a_bar, mat_a_hat, vec_b_b
sum(mat_mu[i, :]) <= vec_b_bar[i])
@constraint(model, vec_lambda[i] .+ hcat(mat_mu[i, :]) .>= hcat(mat_a_hat[i, :]) .* vec_z[:])
end
solve(model)
obj_result = getobjectivevalue(model)
vec_result_x = getvalue(vec_x)
optimize!(model)
obj_result = objective_value(model)
vec_result_x = value(vec_x)
# println("Result: vec_x = $vec_result_x.")
vec_result_z = getvalue(vec_z)
mat_result_mu = getvalue(mat_mu)
vec_result_lambda = getvalue(vec_lambda)
vec_result_z = value(vec_z)
mat_result_mu = value(mat_mu)
vec_result_lambda = value(vec_lambda)
println("---------------------------------------------------------\n",
" 2/2. Nominal Ending\n",
"---------------------------------------------------------\n")
Expand Down
Loading

0 comments on commit e17952e

Please sign in to comment.