Does not converge optimal gas flow with complementarity constraints

[Robbybp]

I'm opening this issue per @amontoison's request.

Here's the matgas file:

gaslib40.m

function mgc = gaslib40

%% required global data
mgc.gas_specific_gravity         = 0.6     ; % dimensionless
mgc.specific_heat_capacity_ratio = 1.4     ; % dimensionless
mgc.temperature                  = 288.706 ; % K
mgc.compressibility_factor       = 1       ; % dimensionless
mgc.units                        = 'si'    ;
%% optional global data
mgc.is_per_unit        = 0     ;
mgc.base_length        = 1000  ; % m
mgc.base_pressure      = 1.0e6 ; % Pa
mgc.base_flow          = 100   ; % kg/s
mgc.economic_weighting = 0.1   ; % dimensionless

%% junction data
% id p_min p_max p_nominal junction_type status pipeline_name edi_id lat lon
mgc.junction = [
32 101325.0 8.101325e6 4.101325e6 0 1 '32.0' '32.0' 7.151456801 47.70262448
29 101325.0 8.101325e6 4.101325e6 0 1 '29.0' '29.0' 6.405925721 48.87915021
1  101325.0 8.101325e6 4.101325e6 0 1 '1.0'  '1.0'  8.462165532 48.68116719
24 101325.0 8.101325e6 4.101325e6 0 1 '24.0' '24.0' 6.985401049 48.99650392
12 101325.0 8.101325e6 4.101325e6 0 1 '12.0' '12.0' 7.444375531 46.93619913
20 101325.0 8.101325e6 4.101325e6 0 1 '20.0' '20.0' 6.732481162 48.45970144
4  101325.0 8.101325e6 4.101325e6 0 1 '4.0'  '4.0'  8.941170283 49.7403095 
2  101325.0 8.101325e6 4.101325e6 0 1 '2.0'  '2.0'  8.583339554 49.33348585
25 101325.0 8.101325e6 4.101325e6 0 1 '25.0' '25.0' 7.375377697 47.06334966
6  101325.0 8.101325e6 4.101325e6 0 1 '6.0'  '6.0'  7.175358607 48.44886966
23 101325.0 8.101325e6 4.101325e6 0 1 '23.0' '23.0' 6.704129879 47.04609847
22 101325.0 8.101325e6 4.101325e6 0 1 '22.0' '22.0' 7.391411724 47.03802643
11 101325.0 8.101325e6 4.101325e6 0 1 '11.0' '11.0' 8.421237101 48.69715728
35 101325.0 8.101325e6 4.101325e6 0 1 '35.0' '35.0' 7.060507549 47.33697273
13 101325.0 8.101325e6 4.101325e6 0 1 '13.0' '13.0' 7.445707841 48.52155545
27 101325.0 8.101325e6 4.101325e6 0 1 '27.0' '27.0' 6.682728687 47.71858108
31 101325.0 8.101325e6 4.101325e6 0 1 '31.0' '31.0' 6.970774678 48.88544272
5  101325.0 8.101325e6 4.101325e6 0 1 '5.0'  '5.0'  8.901176918 49.75682733
15 101325.0 8.101325e6 4.101325e6 0 1 '15.0' '15.0' 6.890115373 48.50659413
33 101325.0 8.101325e6 4.101325e6 0 1 '33.0' '33.0' 7.160019633 47.44368652
28 101325.0 8.101325e6 4.101325e6 0 1 '28.0' '28.0' 8.522974689 49.59906651
16 101325.0 8.101325e6 4.101325e6 0 1 '16.0' '16.0' 9.434658491 48.43337494
14 101325.0 8.101325e6 4.101325e6 0 1 '14.0' '14.0' 7.531026915 48.54893883
39 101325.0 8.101325e6 4.101325e6 0 1 '39.0' '39.0' 6.037965632 48.81954737
40 101325.0 8.101325e6 4.101325e6 0 1 '40.0' '40.0' 8.972793056 49.76190172
38 101325.0 8.101325e6 5.0e6      1 1 '38.0' '38.0' 6.837566076 48.9636071 
21 101325.0 8.101325e6 4.101325e6 0 1 '21.0' '21.0' 7.297282518 47.85672269
7  101325.0 8.101325e6 4.101325e6 0 1 '7.0'  '7.0'  6.076634313 48.81046104
34 101325.0 8.101325e6 4.101325e6 0 1 '34.0' '34.0' 7.012162287 47.32750791
8  101325.0 8.101325e6 4.101325e6 0 1 '8.0'  '8.0'  6.991006806 48.85701639
36 101325.0 8.101325e6 4.101325e6 0 1 '36.0' '36.0' 7.334321825 47.68491579
26 101325.0 8.101325e6 4.101325e6 0 1 '26.0' '26.0' 7.125950441 48.43944074
17 101325.0 8.101325e6 4.101325e6 0 1 '17.0' '17.0' 7.229303087 47.62325261
10 101325.0 8.101325e6 4.101325e6 0 1 '10.0' '10.0' 8.226777359 48.79746338
19 101325.0 8.101325e6 4.101325e6 0 1 '19.0' '19.0' 8.568217723 49.36328555
37 101325.0 8.101325e6 4.101325e6 0 1 '37.0' '37.0' 6.617232238 47.57618394
9  101325.0 8.101325e6 4.101325e6 0 1 '9.0'  '9.0'  6.556714558 46.91925127
18 101325.0 8.101325e6 4.101325e6 0 1 '18.0' '18.0' 6.646695907 47.48827194
30 101325.0 8.101325e6 4.101325e6 0 1 '30.0' '30.0' 6.482711455 48.68585213
3  101325.0 8.101325e6 4.101325e6 0 1 '3.0'  '3.0'  8.5225982   49.35375393
];

%% pipe data
% id fr_junction to_junction diameter length friction_factor p_min p_max status
mgc.pipe = [
32 10 2  0.8 65057.17427 0.002070376 0 1.0e8 1
29 28 5  1.0 32449.37205 0.001988785 0 1.0e8 1
1  38 31 1.0 13071.08523 0.001988785 0 1.0e8 1
24 20 15 1.0 12766.70341 0.001988785 0 1.0e8 1
12 1  16 0.8 76893.55076 0.002070376 0 1.0e8 1
20 15 30 0.8 36061.00987 0.002070376 0 1.0e8 1
4  26 8  0.8 47488.28385 0.002070376 0 1.0e8 1
2  32 21 0.6 20322.20543 0.002183191 0 1.0e8 1
25 30 7  0.8 32921.25982 0.002070376 0 1.0e8 1
6  34 23 0.8 39036.04181 0.002070376 0 1.0e8 1
23 6  13 1.0 21557.56619 0.001988785 0 1.0e8 1
22 30 20 0.8 31179.6191  0.002070376 0 1.0e8 1
11 35 33 0.4 14043.11354 0.002358537 0 1.0e8 1
35 14 10 0.8 58218.96955 0.002070376 0 1.0e8 1
13 33 17 0.6 20634.69827 0.002183191 0 1.0e8 1
27 19 3  0.8 3479.454666 0.001621777 0 1.0e8 1
5  34 35 0.6 3802.586677 0.002183191 0 1.0e8 1
31 10 11 1.0 18136.59729 0.001988785 0 1.0e8 1
15 17 36 0.6 10452.03118 0.002183191 0 1.0e8 1
33 10 3  0.8 65532.21271 0.002070376 0 1.0e8 1
28 4  5  1.0 3418.008251 0.001565017 0 1.0e8 1
16 31 24 0.8 12397.35216 0.002070376 0 1.0e8 1
14 17 32 0.6 10586.12947 0.002183191 0 1.0e8 1
39 27 32 0.6 35218.83909 0.002183191 0 1.0e8 1
21 30 29 0.8 22224.15325 0.002070376 0 1.0e8 1
38 37 18 0.6 10022.78298 0.002183191 0 1.0e8 1
7  35 25 0.4 38659.82436 0.002358537 0 1.0e8 1
34 13 14 1.0 6998.053779 0.001988785 0 1.0e8 1
8  23 9  0.6 18017.8496  0.002183191 0 1.0e8 1
36 26 27 0.8 86690.26557 0.002070376 0 1.0e8 1
26 4  19 0.8 49866.14839 0.002070376 0 1.0e8 1
17 36 21 0.6 19303.19202 0.002183191 0 1.0e8 1
10 22 12 0.4 12015.87483 0.002358537 0 1.0e8 1
19 26 15 1.0 18969.41271 0.001988785 0 1.0e8 1
37 27 37 0.6 16579.326   0.002183191 0 1.0e8 1
9  25 22 0.6 3067.54744  0.002183191 0 1.0e8 1
18 26 21 0.6 66036.59463 0.002183191 0 1.0e8 1
30 28 19 0.8 26427.48165 0.002070376 0 1.0e8 1
3  18 34 0.8 32868.20253 0.002070376 0 1.0e8 1
];

%% compressor data
% id fr_junction to_junction c_ratio_min c_ratio_max power_max flow_min flow_max inlet_p_min inlet_p_max outlet_p_min outlet_p_max status operating_cost
mgc.compressor = [
4 40 4  1 1.5 100000 0 1500 0 1.0e8 0 1.0e8 1 1
1 6  26 1 1.5 100000 0 1500 0 1.0e8 0 1.0e8 1 1
5 39 7  1 1.5 100000 0 1500 0 1.0e8 0 1.0e8 1 1
2 11 1  1 1.5 100000 0 1500 0 1.0e8 0 1.0e8 1 1
6 31 8  1 1.5 100000 0 1500 0 1.0e8 0 1.0e8 1 1
3 19 2  1 1.5 100000 0 1500 0 1.0e8 0 1.0e8 1 1
];

%% receipt data
% id junction_id injection_min injection_max injection_nominal is_dispatchable status offer_price
mgc.receipt = [
1  32 0 0           0 1 1 1   
2  29 0 0           0 1 1 1   
3  24 0 0           0 1 1 1   
4  12 0 0           0 1 1 1   
5  20 0 0           0 1 1 1   
6  25 0 0           0 1 1 1   
7  23 0 0           0 1 1 1   
8  22 0 0           0 1 1 1   
9  11 0 0           0 1 1 1   
10 35 0 0           0 1 1 1   
11 13 0 0           0 1 1 1   
12 27 0 0           0 1 1 1   
13 31 0 0           0 1 1 1   
14 15 0 0           0 1 1 1   
15 33 0 0           0 1 1 1   
16 28 0 0           0 1 1 1   
17 16 0 0           0 1 1 1   
18 14 0 0           0 1 1 1   
19 39 0 8.0902778   0 1 1 1   
20 40 0 158.0902778 0 1 1 1   
21 21 0 0           0 1 1 1   
22 34 0 0           0 1 1 1   
23 36 0 0           0 1 1 1   
24 26 0 0           0 1 1 1   
25 17 0 0           0 1 1 1   
26 10 0 0           0 1 1 1   
27 19 0 0           0 1 1 1   
28 37 0 0           0 1 1 1   
29 9  0 0           0 1 1 1   
30 18 0 0           0 1 1 1   
31 30 0 0           0 1 1 1   
32 38 0 1300        0 1 1 0.75
];

%% delivery data
% id junction_id withdrawal_min withdrawal_max withdrawal_nominal is_dispatchable status bid_price
mgc.delivery = [
1  32 0 16.35416667 0 1 1 3   
2  29 0 16.35416667 0 1 1 3   
3  24 0 16.35416667 0 1 1 3   
4  12 0 16.35416667 0 1 1 3   
5  20 0 36.35416667 0 1 1 3   
6  25 0 1.35416667  0 1 1 3   
7  23 0 16.35416667 0 1 1 3   
8  22 0 1.35416667  0 1 1 3   
9  11 0 16.35416667 0 1 1 3   
10 35 0 16.35416667 0 1 1 3   
11 13 0 16.35416667 0 1 1 3   
12 27 0 46.35416667 0 1 1 3   
13 31 0 16.35416667 0 1 1 3   
14 15 0 16.35416667 0 1 1 3   
15 33 0 36.35416667 0 1 1 3   
16 28 0 16.35416667 0 1 1 3   
17 16 0 36.35416667 0 1 1 3   
18 14 0 16.35416667 0 1 1 3   
19 39 0 0           0 1 1 3   
20 40 0 0           0 1 1 3   
21 21 0 24.35416667 0 1 1 3   
22 34 0 16.35416667 0 1 1 3   
23 36 0 16.35416667 0 1 1 3   
24 26 0 22.35416667 0 1 1 3   
25 17 0 16.35416667 0 1 1 3   
26 10 0 0           0 1 1 3   
27 19 0 45.35416667 0 1 1 3   
28 37 0 16.35416667 0 1 1 3   
29 9  0 43.35416667 0 1 1 3   
30 18 0 16.35416667 0 1 1 3   
31 30 0 105.3541667 0 1 1 3   
32 38 0 0           0 1 1 2.75
];

end

Here's the code:

import GasModels
import NLPModelsJuMP
import MadNCL
import MadNLPHSL
fpath = "gaslib40.m"
data = GasModels.parse_file(fpath)
gm = GasModels.instantiate_model(data, GasModels.CWPGasModel, GasModels.build_ogf)
nlp = NLPModelsJuMP.MathOptNLPModel(gm.model)
res = MadNCL.madncl(nlp; linear_solver = MadNLPHSL.Ma57Solver)

FWIW, I can solve this problem with a non-differentiable formulation (GasModels.WPGasModel) using SLP. I have not yet verified that this solution is also a solution for the complementarity formulation (although they should be equivalent).

[frapac]

Thank you for posting a detailed MWE. I confirm I am able to reproduce the issue on my machine.

It turns out that this nonlinear program is very sensitive to the linear solver we are using (HSL ma27 and HSL ma57 are returning different results). Apparently the KKT system K2rAuglagKKTSystem gives better results than the default KKT system. I am able to solve the problem to tol=1e-5 using the following options:

fpath = "scripts/gaslib40.m"
data = GasModels.parse_file(fpath)
gm = GasModels.instantiate_model(data, GasModels.CWPGasModel, GasModels.build_ogf)
nlp = NLPModelsJuMP.MathOptNLPModel(gm.model)

ncl_options = MadNCL.NCLOptions{Float64}(;
    opt_tol=1e-5,   # relax tolerance
    feas_tol=1e-5,   # relax tolerance
    scaling_max_gradient=100.0,  # use Ipopt's scaling
)
res = @time MadNCL.madncl(
    nlp;
    ncl_options=ncl_options,
    linear_solver=MadNLPHSL.Ma57Solver,
    kkt_system=MadNCL.K2rAuglagKKTSystem,
    print_level=MadNLP.ERROR,
)

I obtain the following convergence:

MadNCL algorithm

Total number of variables............................:      240
Total number of constraints..........................:      285

outer  inner     objective    inf_pr   inf_du    η        μ       ρ 
    0      0 +0.0000000e+00 0.00e+00 0.00e+00 1.00e-01 1.0e-01 1.00e+02
    1r    12 -5.1042843e+03 7.82e-01 3.59e+05 1.00e-02 1.0e-01 1.00e+03
    2r    14 -5.0954960e+03 7.85e-01 3.58e+05 1.00e-02 1.0e-01 1.00e+04
    3r    27 -5.0458290e+03 7.92e-01 2.72e+07 1.00e-02 1.0e-01 1.00e+05
    4r    29 -5.0446682e+03 7.92e-01 3.81e+07 1.00e-02 1.0e-01 1.00e+06
    5r    37 -4.9212159e+03 7.72e-01 5.60e+07 1.00e-02 1.0e-01 1.00e+07
    6r    38 -4.9212159e+03 7.72e-01 5.60e+07 1.00e-02 1.0e-01 1.00e+08
    7r    39 -4.9212159e+03 7.72e-01 7.72e+07 1.00e-02 1.0e-01 1.00e+09
    8    109 -1.4432182e+02 1.92e-05 2.40e+01 1.02e-03 1.0e-02 1.00e+09
    9    112 -1.4844335e+02 1.91e-06 6.55e+01 1.10e-05 1.1e-04 1.00e+09
   10    140 -1.4880998e+02 1.06e-05 6.41e-04 1.32e-09 1.3e-08 1.00e+09
   11    338 -1.4881412e+02 5.57e-07 8.37e-06 1.00e-10 1.0e-09 1.00e+09
  0.077764 seconds (289.36 k allocations: 28.022 MiB)
Optimal solution found.

@amontoison we might want to change the option to use K2rAuglagKKTSystem by default.

@Robbybp Your instance is interesting. You implement the complementarity constraints using the following product: f_plus * f_minus = 0.0. This can lead to some numerical issue in the linear solver. In general we recommend rewriting the complementarity constraint as f_plus * f_minus <= 0 if you solve the problem with IPM.

Now that the nonlinear functions are first-class citizens in JuMP, we are currently developing a new MOI wrapper to automatically reformulate the complementarity constraint: https://github.com/blegat/ComplementOpt.jl
It would be interesting to rewrite this instance in GasModels.jl with MOI.Complements, and test the solution we obtain using the different reformulations implemented in GasModels.jl.

If you are interested, we are currently gathering a collection of large-scale MPCCs in this repo: https://github.com/MadNLP/MPCCbenchmark.jl . We gather instances from optimal control (with switch detection), from equilibrium problems (MacEPEC) and from power systems (SCOPF, power flow with PV/PQ switches). If you have other problems with complementarity constraints, I would be interested to add them in this new benchmark. Let me know if you are interested!

[amontoison]

@frapac Is K2 with the SparseKKTSystem formulation is the default one?
Yes, we should define K2r as default, it is why we did MadNCL! :)