demand function is not working

[floswald]

this call is outdated

https://github.com/ScPo-CompEcon/CourseMatch.jl/blob/master/src/demand.jl#L117

update to new API. you need to do ] up in teh coursematch project

[jeannesorin]

Does m = Model(with_optimizer(Ipopt.Optimizer))work or is there an alternative for CbcSolver? I could not find anything else online.

[floswald]

   _       _ _(_)_     |  Documentation: https://docs.julialang.org
  (_)     | (_) (_)    |
   _ _   _| |_  __ _   |  Type "?" for help, "]?" for Pkg help.
  | | | | | | |/ _` |  |
  | | |_| | | | (_| |  |  Version 1.1.0 (2019-01-21)
 _/ |\__'_|_|_|\__'_|  |  Official https://julialang.org/ release
|__/                   |

julia> using JuMP

julia> using Cbc

julia> m = Model(with_optimizer(Cbc.Optimizer))
A JuMP Model
Feasibility problem with:
Variables: 0
Model mode: AUTOMATIC
CachingOptimizer state: EMPTY_OPTIMIZER
Solver name: COIN Branch-and-Cut (Cbc)

[floswald]

anyone working on this?

[danielbgyetvai]

Yes, I am. I have updated the constraint calls in demand!() and am now trying to figure out how to adapt the objective.

[floswald]

excellent! keep us posted with how it goes / any questions.

[danielbgyetvai]

It is not directly relevant but I created the following function (replacing and adapting rand_constraints() which was not working) in student.jl:

#This is an auxiliary function for testing purposes only. It might create inconsistencies, economically speaking.
function add_rand_constraints!(s::Student)
    N = length(s.constraints.isFC)
    if N<2
        error("N must be > 2")
    end
    s.constraints.mandatory = rand([true false],N)
    s.constraints.notTCprogram = rand([true false false false false false],N)
    s.constraints.notTCsemester = rand([true false false false false false],N)
    s.constraints.isTC = rand([true false],N)
    s.constraints.notFCprogram = rand([true false false false false false],N)
    s.constraints.notFCsemester = rand([true false false false false false],N)
    s.constraints.isFC = rand([true false],N)
    s.constraints.notELprogram = rand([true false false false false false],N)
    s.constraints.notELsemester = rand([true false false false false false],N)
    s.constraints.isEL = rand([true false],N)
end 

This can still be improved since I feel it might create economically inconsistent sets of constraints like courses for which isTC is 0 and notTCprogram is 0 as well.

[gdannay]

Hey floswald,

We have some troubles with the solver. We manage to change all the constraints but the Cbc solver does not handle NL objective function. We tried with Ipopt but it doesn't handle mixed integers. Do you know a solver that can solve Mixed Integers Non Linear Programming problems ?

[floswald]

do we have a nonlinear objective?

[gdannay]

The objective is @objective(m, Max, x's.preferences.x)

[floswald]

Yes. That’s a linear function. Not?

On Tue 30 Apr 2019 at 17:15, Grégory Dannay notifications@github.com wrote:

The objective is @objective https://github.com/objective(m, Max, x's.preferencesx)

— You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub https://github.com/ScPo-CompEcon/CourseMatch.jl/issues/1#issuecomment-487993550, or mute the thread https://github.com/notifications/unsubscribe-auth/AAHZM5UNU4KND62BVXQJYI3PTBPCZANCNFSM4HEGY4SQ .

[floswald]

i think i know where you are coming from though, so this is a good question. of course you are right that in theory this defines a quadratic form, so non-linear. notice that if we don't multiply it out, it's linear. what may be even more helpful though is to go back to what your predecessors did. x' * s.preferences * x is the assumed utility function of a student. the whole purpose of this is to be able to deal with a matrix of preferences. this is necessary if you want to allow to define complementarities. I really want course A but only if I can also get course B at the same time. This would be indicated by a high value in cell (A,B) of matrix s.preferences. note finally that the quadratic form is really just the sum of utilities (given x \in {0,1}!).

[danielbgyetvai]

Ok, so our task is to transform the objective such that it becomes linear in the elements of x. (And not to look for an optimizer that would solve the problem in its current form.) By the way, with @gdannay and @SimonGuillard we noticed that the clash constraint requires a similar treatment since that one is also quadratic in x.

[floswald]

well, that formulation used to work, so I am not sure what the problem is with this. are you saying that you get an error message if you use that form of objective? (again, feel free to assume another utility function, but be careful with what you choose here - simple sum of utilities seems like a solid choice!)

[danielbgyetvai]

I am not sure it would work given the abs() or that it is the simplest way but I have found the following: if x[i] and x[j] are binary (0 or 1) then x[i]*x[j] = ((x[i] + x[j]) - abs(x[i] - x[j]))/2.

[floswald]

that's not what i meant. I meant you should just see if this prototype works with Cbc instead of Gurobi (and the up do date JuMP API)

using JuMP, Gurobi

m = Model(solver = GurobiSolver())

@variable(m, x[1:3], Bin)

utility = [80 10 -50; 10 90 0; -50 0 60]
price = [20, 30, 15]
budget = 50

@objective(m, Max, x'*utility*x)

#profit = [10 20 30]

@constraint(m, price'*x <= budget)

status = solve(m)

for i in 1:3
    println("x[$i] = ", getvalue(x[i]))
end

if it does, then - why change it?

[floswald]

indeed the quadratic form doesn't work with Cbc. We can fall back to gurobi. this works for me in the meantime:

julia> using JuMP, Cbc
[ Info: Recompiling stale cache file /Users/florian.oswald/.julia/compiled/v1.1/JuMP/DmXqY.ji for JuMP [4076af6c-e467-56ae-b986-b466b2749572]
[ Info: Recompiling stale cache file /Users/florian.oswald/.julia/compiled/v1.1/Cbc/ARPfV.ji for Cbc [9961bab8-2fa3-5c5a-9d89-47fab24efd76]

julia> m = Model(with_optimizer(Cbc.Optimizer))
A JuMP Model
Feasibility problem with:
Variables: 0
Model mode: AUTOMATIC
CachingOptimizer state: EMPTY_OPTIMIZER
Solver name: COIN Branch-and-Cut (Cbc)

julia> @variable(m, x[1:3], Bin)
3-element Array{VariableRef,1}:
 x[1]
 x[2]
 x[3]

julia> utility = [80 10 -50; 10 90 0; -50 0 60]
3×3 Array{Int64,2}:
  80  10  -50
  10  90    0
 -50   0   60

julia> price = [20, 30, 15]
3-element Array{Int64,1}:
 20
 30
 15

julia> budget = 50
50

julia> exp = @expression(m, sum(x[i] * utility[i] for i in 1:3))
80 x[1] + 10 x[2] - 50 x[3]

julia> @objective(m,Max, exp)
80 x[1] + 10 x[2] - 50 x[3]

julia> @constraint(m,sum(x[i] * price[i] for i in 1:3) <= budget)
20 x[1] + 30 x[2] + 15 x[3] ≤ 50.0

julia> optimize!(m)
Welcome to the CBC MILP Solver 
Version: 2.9.9 
Build Date: Dec 31 2018 

command line - Cbc_C_Interface -solve -quit (default strategy 1)
Continuous objective value is 90 - 0.00 seconds
Cgl0004I processed model has 0 rows, 0 columns (0 integer (0 of which binary)) and 0 elements
Cbc3007W No integer variables - nothing to do
Cuts at root node changed objective from -90 to -1.79769e+308
Probing was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds)
Gomory was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds)
Knapsack was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds)
Clique was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds)
MixedIntegerRounding2 was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds)
FlowCover was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds)
TwoMirCuts was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds)

Result - Optimal solution found

Objective value:                90.00000000
Enumerated nodes:               0
Total iterations:               0
Time (CPU seconds):             0.00
Time (Wallclock seconds):       0.01

Total time (CPU seconds):       0.00   (Wallclock seconds):       0.02

julia> for i in 1:3
       println("x[$i] = $(value(x[i]))")
       end
x[1] = 1.0
x[2] = 1.0
x[3] = 0.0

[floswald]

installation of gurobi only works while you are on the ScPo wifi. follow here for gurobi https://github.com/JuliaOpt/Gurobi.jl - you need to get a free license at gurobi.com (as said, for first activation need to be at ScPo). your decision if we need the quadratic form, I think it's a nice way to deal with the complementary preferences.

[floswald]

Yeah I remember now that this quadratic obj func was the reason we went with gurobi. I guess the benefits of this simple representation of cross preferences is a big advantage so I would vote to keep gurobi. Thoughts?

[floswald]

alright, so Gurobi works. given there were no objections to that, let's go back to Gurobi

julia> using JuMP, Gurobi

julia> m = Model(with_optimizer(Gurobi.Optimizer))
Academic license - for non-commercial use only
A JuMP Model
Feasibility problem with:
Variables: 0
Model mode: AUTOMATIC
CachingOptimizer state: EMPTY_OPTIMIZER
Solver name: Gurobi

julia> @variable(m, x[1:3], Bin)
3-element Array{VariableRef,1}:
 x[1]
 x[2]
 x[3]

julia> utility = [80 10 -50; 10 90 0; -50 0 60]
3×3 Array{Int64,2}:
  80  10  -50
  10  90    0
 -50   0   60

julia> price = [20, 30, 15]
3-element Array{Int64,1}:
 20
 30
 15

julia> budget = 50
50

julia> @objective(m, Max, x'*utility*x)
80 x[1]² + 20 x[1]*x[2] - 100 x[1]*x[3] + 90 x[2]² + 60 x[3]²

julia> #profit = [10 20 30]

       @constraint(m, price'*x <= budget)
20 x[1] + 30 x[2] + 15 x[3] ≤ 50.0

julia> optimize!(m)
Academic license - for non-commercial use only
Optimize a model with 1 rows, 3 columns and 3 nonzeros
Model has 5 quadratic objective terms
Variable types: 0 continuous, 3 integer (3 binary)
Coefficient statistics:
  Matrix range     [2e+01, 3e+01]
  Objective range  [0e+00, 0e+00]
  QObjective range [4e+01, 2e+02]
  Bounds range     [0e+00, 0e+00]
  RHS range        [5e+01, 5e+01]
Found heuristic solution: objective -0.0000000
Presolve time: 0.00s
Presolved: 4 rows, 5 columns, 10 nonzeros
Variable types: 0 continuous, 5 integer (5 binary)

Root relaxation: objective 1.900000e+02, 0 iterations, 0.00 seconds

    Nodes    |    Current Node    |     Objective Bounds      |     Work
 Expl Unexpl |  Obj  Depth IntInf | Incumbent    BestBd   Gap | It/Node Time

*    0     0               0     190.0000000  190.00000  0.00%     -    0s

Explored 0 nodes (0 simplex iterations) in 0.00 seconds
Thread count was 4 (of 4 available processors)

Solution count 2: 190 -0 

Optimal solution found (tolerance 1.00e-04)
Best objective 1.900000000000e+02, best bound 1.900000000000e+02, gap 0.0000%

julia> for i in 1:3
           println("x[$i] = $(value(x[i]))")
       end
x[1] = 1.0
x[2] = 0.9999999999999958
x[3] = 0.0