How to model an implication

Specific Domains Optimization (Mathematical)
21

This very simplified form also gives similar results.

julia>  main() do n, N, model, x
                   @variable(model, z, Int)
                   @constraint(model, z == sum(x[1, :]) - N)   
                   @constraint(model, !z --> {sum(x[2, :]) == 2})
               end
Test Summary: | Pass  Total  Time
main          | 4070   4070  0.8s

Is it just a coincidence? Or are the formulations(*) actually equivalent?

(*)I’m having trouble following the logic of the second script.
For example, regarding the statement
@constraint(model, z[3] == z[1] + z[2])

I understand that if z[3]==0 then z[1]==z[2]==0.

I don’t understand how it works in different cases.
For example, if it makes sense, since z is Bin, what happens in the case z[1]==z[2]==1?

22
Clarabel.jl Apache LP, QP, SOCP, SDP
COSMO.jl Apache LP, QP, SOCP, SDP
DAQP DAQP.jl MIT (Mixed-binary) QP
HiGHS HiGHS.jl MIT (MI)LP, QP
Ipopt Ipopt.jl EPL LP, QP, NLP
MadNLP.jl MIT LP, QP, NLP
NLopt NLopt.jl GPL LP, QP, NLP
OSQP OSQP.jl Apache LP, QP
SCS SCS.jl MIT LP, QP, SOCP, SDP
StatusSwitchingQP.jl MIT LP, QP

These are the ones that deal with quadratic problems.
Which one could I use for my model, with a quadratic objective function?

23

Is it just a coincidence?

I need to take a better look at this. I would have expected an error. The z variable in an indicator must be binary.

For example, if it makes sense, since z is Bin, what happens in the case z[1]==z[2]==1?

This can’t happen. The problem would be infeasible. In other words, you can’t have both s[1] and s[2] be non-zero at the same time.

Which one could I use for my model, with a quadratic objective function?

If you want a quadratic objective, you should choose SCIP or Gurobi (or maybe DAQP?).

Alternatively, if all the variables in the quadratic, you can linearise the terms as follows:

model = Model()
@variable(model, x[1:2], Bin)
# @objective(model, Max, x[1] * x[2]) ?
@variable(model, z, Bin)
@constraints(model, begin
    z <= x[1]
    z <= x[2]
    z >= x[1] + x[2] - 1
end)
@objective(model, Max, z)
24

If this helps your analysis, consider the following fact (*).


model = Model(HiGHS.Optimizer)
# set_silent(model)
@variable(model, x[1:nr, 1:ng, turni ], Bin)
@variable(model,0<=s[1:nr]<=nr, Int)
@constraints(model, begin

    # each day must have exactly one Ti or one Fi
    fer1[k in union(t_feriali,t_card_urg), g in g_feriali], sum(x[:,g,k]) == 1
    fer2[k in t_festivi, g in g_feriali], sum(x[:,g,k]) == 0

    fest1[k in t_festivi, g in g_festivi], sum(x[:,g,k]) == 1
    fest2[k in union(t_feriali, t_card_urg), g in g_festivi], sum(x[:,g,k]) == 0
    
    # range of N and NF for row for month
    [r in ops_lim], sum(x[r,j,"N"]+x[r,j,"NF"] for j in 1:ng) <= 4
    [r in ops_lim], sum(x[r,j,"N"]+x[r,j,"NF"] for j in 1:ng) >= 3
    
    [r in ops_lim], sum(x[r,j,"N"] for j in 1:ng) <= 3
    [r in ops_lim], sum(x[r,j,"N"] for j in 1:ng) >= 1
   
    # # NF==2 ==> N<=1
 
      [r in ops_lim],  s[r] == sum(x[r,:,"NF"])-2
      [r in ops_lim], !s[r] --> {sum(x[r,:,"N"]) <= 1}
   ...
# many other constraints
end

...
solution_summary(; result = 1, verbose = false)
β”œ solver_name          : HiGHS
β”œ Termination
β”‚ β”œ termination_status : INFEASIBLE
β”‚ β”œ result_count       : 1
β”‚ β”œ raw_status         : kHighsModelStatusInfeasible
β”‚ β”” objective_bound    : 0.00000e+00
β”œ Solution (result = 1)
β”‚ β”œ primal_status        : NO_SOLUTION
β”‚ β”œ dual_status          : INFEASIBILITY_CERTIFICATE
β”‚ β”œ objective_value      : 0.00000e+00
β”‚ β”œ dual_objective_value : 6.00000e+00
β”‚ β”” relative_gap         : Inf

model = Model(HiGHS.Optimizer)
# set_silent(model)
@variable(model, x[1:nr, 1:ng, turni ], Bin)
@variable(model, z[1:nr, 1:ng], Bin)
@constraints(model, begin

    # each day must have exactly one Ti or one Fi
    fer1[k in union(t_feriali,t_card_urg), g in g_feriali], sum(x[:,g,k]) == 1
    fer2[k in t_festivi, g in g_feriali], sum(x[:,g,k]) == 0

    fest1[k in t_festivi, g in g_festivi], sum(x[:,g,k]) == 1
    fest2[k in union(t_feriali, t_card_urg), g in g_festivi], sum(x[:,g,k]) == 0
  
    
     # range of N and NF for row for month
    [r in ops_lim], sum(x[r,j,"N"]+x[r,j,"NF"] for j in 1:ng) <= 4
    [r in ops_lim], sum(x[r,j,"N"]+x[r,j,"NF"] for j in 1:ng) >= 3
    
    [r in ops_lim], sum(x[r,j,"N"] for j in 1:ng) <= 3
    [r in ops_lim], sum(x[r,j,"N"] for j in 1:ng) >= 1
   
    # # NF==2 ==> N<=1
      [r in ops_lim], sum(x[r,:,"NF"]) == sum(i * z[r, i] for i in 1:ng)
      [r in ops_lim], z[r,2] --> {sum(x[r, :, "N"]) <= 1}
      [r in ops_lim], sum(z[r,:]) <= 1

...

Solving report
  Status            Optimal
  Primal bound      0
  Dual bound        0
  Gap               0% (tolerance: 0.01%)
  P-D integral      0
  Solution status   feasible        
                    0 (objective)   
                    0 (bound viol.) 
                    0 (int. viol.)  
                    0 (row viol.)

(*)
I tried to apply one of your suggestions to my case, but I’m not entirely sure if I did it right

PS

The following form also passes the parser checks and provides results consistent with the intended purpose of the constraint

    # # NF==2 ==> N<=1 for r in ops_lim
 
    s[ops_lim] .== sum(x[ops_lim,:,"NF"])-2

    [r in ops_lim], !s[r] --> {sum(x[r,:,"N"]) <= 1}

The following form, however, gives an unfeasible result

# NF==2 ==> N<=1 for r in ops_lim
 
   [r in ops_lim],   s[r] == sum(x[r,:,"NF"])-2

    [r in ops_lim], !s[r] --> {sum(x[r,:,"N"]) <= 1}
25

Indicator constraints are not valid when the indicator variable is not binary. We should be throwing an error here: [Bridges] no error thrown when IndicatorToMILPBridge does not imply first variable is binary Β· Issue #2856 Β· jump-dev/MathOptInterface.jl Β· GitHub