Problems with a picewise quadratic problem


#1

Hi Guys,

I need to model the following optimization problem

\sum_{j=1}^{N} (\alpha * w_j (d_j(x)-p_j)^2)

where \alpha = 1 when (d_j(x)-p_j) > 0 and \alpha = 0 otherwise

I’ve been through the JuMP’s documentation without much success. Could anyone tell me whether JuMP supports this kind of mathematical expressions?

Thanks!


#2

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 https://github.com/joehuchette/PiecewiseLinearOpt.jl which will probably give you a more efficient formulation, but I’m not as familiar with its interface.


#3

Thank you very much @rdeits!!

Cheers

Guillermo