ConditionalJuMP.jl is a tool that I wrote to make things like your problem easier to formulate in JuMP. For example, if we introduce a new variable z
and constrain that z_j = d_j - p_j
when d_j - p_j >= 0
and z_j = 0
otherwise, then we have the same problem you wrote, but with z = alpha * (d_j - p_j)
.
using JuMP, ConditionalJuMP, Gurobi
N = 10
w = randn(N)
p = randn(N)
model = Model(solver=GurobiSolver())
@variable model d[1:N] lowerbound=-10 upperbound=10
@variable model z[1:N] lowerbound=0 upperbound=100
for i in 1:N
@implies(model,
(d[i] - p[i] >= 0) => z[i] == d[i] - p[i],
(d[i] - p[i] <= 0) => z[i] == 0)
end
@objective model Min sum(w[i] * z[i]^2 for i in 1:N)
solve(model)
@show getvalue(d) getvalue(z);
gives:
Optimize a model with 80 rows, 30 columns and 190 nonzeros
Model has 10 quadratic objective terms
Variable types: 20 continuous, 10 integer (10 binary)
Coefficient statistics:
Matrix range [1e+00, 1e+02]
Objective range [0e+00, 0e+00]
QObjective range [1e-01, 5e+00]
Bounds range [1e+00, 1e+02]
RHS range [1e+01, 1e+02]
Found heuristic solution: objective 0
Presolve removed 80 rows and 30 columns
Presolve time: 0.00s
Presolve: All rows and columns removed
Explored 0 nodes (0 simplex iterations) in 0.00 seconds
Thread count was 1 (of 12 available processors)
Solution count 1: -732.788
Pool objective bound -732.788
Optimal solution found (tolerance 1.00e-04)
Best objective -7.327877477104e+02, best bound -7.327877477104e+02, gap 0.0000%
getvalue(d) = [10.0, 1.46984, 10.0, 10.0, -2.17514, 0.285866, 10.0, 10.0, -0.522088, 10.0]
getvalue(z) = [10.8604, 0.0, 10.827, 9.46254, 0.0, 0.0, 9.21451, 11.013, 0.0, 9.64365]
or, equivalently:
model = Model(solver=GurobiSolver())
@variable model d[1:N] lowerbound=-10 upperbound=10
@variable model z[1:N] lowerbound=0 upperbound=100
for i in 1:N
@disjunction(model,
(d[i] - p[i] >= 0) & (z[i] == d[i] - p[i]),
(z[i] == 0))
end
@objective model Min sum(w[i] * z[i]^2 for i in 1:N)
solve(model)
@show getvalue(d) getvalue(z);
Note that ConditionalJuMP requires that all your decision variables have bounds.
There’s also GitHub - joehuchette/PiecewiseLinearOpt.jl: Solve optimization problems containing piecewise linear functions which will probably give you a more efficient formulation, but I’m not as familiar with its interface.