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


#4

Hi @rdeits.
It seems that ConditionalJuMP.jl is not supported with Julia 0.7.
Do you plan to fix it, or do you have any advice to do other?
Thanks in advance,
Arnaud.


#5

Someday? Yes, definitely. But probably not in the next month at least. Porting it to 0.7 would really mean porting it to use the new MathOptInterface, which would probably involve rewriting the package to some extent. I think that’s definitely worth doing, but it’s not my top priority right now.


#6

Understood. Thanks Robin for your prompt answer.
I’ll carefully check the availability as I’m a strong user of your awesome package.


#7

Hi guys,

Does anyone knows how to formulate a JuMP Variable comparison?
I was wondering how to fomulate this, just with a if else statement, like:

@variable(model, 0 <= loc_pn_compat_free[loc = LOCATIONS, pn = PN_ALL, time = TIMES] <= 1, Bin)
@variable(model, 0 <= IsInLocked[loc = LOCATIONS, pn = PN_ALL, time = TIMES] <= 1, Bin)

for loc in LOCATIONS, pn in PN_ALL, time in TIMES
	if loc_pn_compat_free[loc,pn,time] >= 0.5
		@constraint(model,IsInLocked[loc,pn,time] == 0)
	end
end

By doing this, I obtain:

ERROR: MethodError: no method matching isless(::JuMP.Variable, ::Float64)

(MILP - Cbc Solver)

The @implies macro does the job great, but only for 1 statement to compare.
Basically, I’d like to simulate something like: @implies(model, x >= 0.5 && y <= 1 => …), which is impossible with ConditionalJuMP package, but maybe possible by another way …?

Thanks in advance,
Arnaud.