[WIP] Experimental SLP algorithm

[Robbybp]

This is a very basic implementation that will need some improvements before production use, but this is probably a good point to get some comments on my design. The algorithm converges several of the Matgas files I've written from Luke and Anatoly's CSV data, but does not converge any of the Matgas files I tried from this repo's tests.

Development roadmap:

• Handle components that did not show up in test data, e.g., valve, regulator, short-pipe, etc.
• Trust region
• Custom initialization options (beyond just using start_value)
• Potentially, a post-optimization Newton solve to clean up primal feasibility

Here's an example that works:

model30ss.m

function mgc = model30ss

%% 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

%% junction data
% id p_min p_max p_nominal junction_type status pipeline_name edi_id lat lon
mgc.junction = [
1  3.447378645e6 5.515805832e6 3.547378645e6  1 1 '1.0'  '1.0'  -0.755  0.8714 
2  3.447378645e6 5.515805832e6 4.4815922385e6 0 1 '2.0'  '2.0'  -0.2421 0.4018 
3  3.447378645e6 5.515805832e6 4.4815922385e6 0 1 '3.0'  '3.0'  0.1248  0.3762 
4  3.447378645e6 5.515805832e6 4.4815922385e6 0 1 '4.0'  '4.0'  0.2794  0.4099 
5  3.447378645e6 5.515805832e6 4.4815922385e6 0 1 '5.0'  '5.0'  0.2794  0.6983 
6  3.447378645e6 5.515805832e6 4.4815922385e6 0 1 '6.0'  '6.0'  0.139   0.8161 
7  3.447378645e6 5.515805832e6 4.4815922385e6 0 1 '7.0'  '7.0'  0.6232  0.72   
8  3.447378645e6 5.515805832e6 4.4815922385e6 0 1 '8.0'  '8.0'  0.7751  0.8401 
9  3.447378645e6 5.515805832e6 4.4815922385e6 0 1 '9.0'  '9.0'  -0.245  0.1949 
10 3.447378645e6 5.515805832e6 4.4815922385e6 0 1 '10.0' '10.0' -0.6576 -0.1382
11 3.447378645e6 5.515805832e6 4.4815922385e6 0 1 '11.0' '11.0' -0.8639 -0.1502
12 3.447378645e6 5.515805832e6 4.4815922385e6 0 1 '12.0' '12.0' -0.9126 0.0901 
13 3.447378645e6 5.515805832e6 4.4815922385e6 0 1 '13.0' '13.0' -0.9069 -0.3377
14 3.447378645e6 5.515805832e6 4.4815922385e6 0 1 '14.0' '14.0' -0.3064 -0.2066
15 3.447378645e6 5.515805832e6 4.4815922385e6 0 1 '15.0' '15.0' -0.2564 -0.4555
16 3.447378645e6 5.515805832e6 4.4815922385e6 0 1 '16.0' '16.0' -0.4456 -0.7103
17 3.447378645e6 5.515805832e6 4.4815922385e6 0 1 '17.0' '17.0' -0.6777 -0.7656
18 3.447378645e6 5.515805832e6 4.4815922385e6 0 1 '18.0' '18.0' -0.7092 -0.869 
19 3.447378645e6 5.515805832e6 4.4815922385e6 0 1 '19.0' '19.0' -0.2994 -0.7993
20 3.447378645e6 5.515805832e6 4.4815922385e6 0 1 '20.0' '20.0' 0.2909  -0.1264
21 3.447378645e6 5.515805832e6 4.4815922385e6 0 1 '21.0' '21.0' 0.3539  -0.3353
22 3.447378645e6 5.515805832e6 4.4815922385e6 0 1 '22.0' '22.0' 0.6117  -0.4363
23 3.447378645e6 5.515805832e6 4.4815922385e6 0 1 '23.0' '23.0' 0.7436  -0.4531
24 3.447378645e6 5.515805832e6 4.4815922385e6 0 1 '24.0' '24.0' 0.7019  -0.9228
25 3.447378645e6 5.515805832e6 4.4815922385e6 0 1 '25.0' '25.0' 0.6146  -0.2656
26 3.447378645e6 5.515805832e6 4.4815922385e6 0 1 '26.0' '26.0' -0.685  0.8014 
27 3.447378645e6 5.515805832e6 4.4815922385e6 0 1 '27.0' '27.0' -0.2421 0.3218 
28 3.447378645e6 5.515805832e6 4.4815922385e6 0 1 '28.0' '28.0' 0.1848  0.3762 
29 3.447378645e6 5.515805832e6 4.4815922385e6 0 1 '29.0' '29.0' -0.3064 -0.2866
30 3.447378645e6 5.515805832e6 4.4815922385e6 0 1 '30.0' '30.0' 0.3009  -0.1764
];

%% pipe data
% id fr_junction to_junction diameter length friction_factor p_min p_max status
mgc.pipe = [
1  26 2  0.9144 100000 0.01 0 1.0e8 1
2  2  3  0.635  30000  0.01 0 1.0e8 1
3  28 4  0.635  5000   0.01 0 1.0e8 1
4  4  5  0.635  15000  0.01 0 1.0e8 1
5  5  6  0.635  10000  0.01 0 1.0e8 1
6  5  7  0.635  5000   0.01 0 1.0e8 1
7  7  8  0.635  10000  0.01 0 1.0e8 1
8  27 9  0.9144 5000   0.01 0 1.0e8 1
9  9  10 0.9144 60000  0.01 0 1.0e8 1
10 10 11 0.635  5000   0.01 0 1.0e8 1
11 11 12 0.635  8000   0.01 0 1.0e8 1
12 11 13 0.635  6000   0.01 0 1.0e8 1
13 10 14 0.9144 80000  0.01 0 1.0e8 1
14 29 15 0.9144 10000  0.01 0 1.0e8 1
15 15 16 0.9144 20000  0.01 0 1.0e8 1
16 16 17 0.635  3000   0.01 0 1.0e8 1
17 17 18 0.635  6000   0.01 0 1.0e8 1
18 16 19 0.635  5000   0.01 0 1.0e8 1
19 15 20 0.9144 40000  0.01 0 1.0e8 1
20 30 21 0.9144 5000   0.01 0 1.0e8 1
21 21 22 0.9144 20000  0.01 0 1.0e8 1
22 22 23 0.9144 5000   0.01 0 1.0e8 1
23 23 24 0.9144 16000  0.01 0 1.0e8 1
24 22 25 0.635  8000   0.01 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 = [
1 1  26 1 1.4 2609950 0 188.2844812 0 1.0e8 0 1.0e8 1 1
2 2  27 1 1.4 1864250 0 164.243841  0 1.0e8 0 1.0e8 1 1
3 3  28 1 1.4 1118550 0 116.1625607 0 1.0e8 0 1.0e8 1 1
4 14 29 1 1.4 745700  0 164.243841  0 1.0e8 0 1.0e8 1 1
5 20 30 1 1.4 745700  0 164.243841  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  6  0 10  0 1 1 12  
2  8  0 26  0 1 1 10  
3  24 0 10  0 1 1 13  
4  25 0 17  0 1 1 12  
5  12 0 22  0 1 1 13  
6  13 0 24  0 1 1 10.5
7  18 0 20  0 1 1 13.5
8  19 0 15  0 1 1 16  
9  12 0 30  0 1 1 6   
10 25 0 0   0 1 1 0   
11 6  0 0   0 1 1 0   
12 18 0 0   0 1 1 0   
13 24 0 0   0 1 1 0   
14 24 0 0   0 1 1 0   
15 25 0 0   0 1 1 0   
16 1  0 110 0 1 1 2   
];

%% delivery data
% id junction_id withdrawal_min withdrawal_max withdrawal_nominal is_dispatchable status bid_price
mgc.delivery = [
1  6  0 0 18   0 1 0   
2  8  0 0 16.5 0 1 0   
3  24 0 0 15   0 1 0   
4  25 0 0 13.5 0 1 0   
5  12 0 0 16.5 0 1 0   
6  13 0 0 13.5 0 1 0   
7  18 0 0 15   0 1 0   
8  19 0 0 18   0 1 0   
9  12 0 0 7.5  0 1 0   
10 25 0 3 6    0 1 20  
11 6  0 4 7    0 1 16  
12 18 0 5 5.5  0 1 13.5
13 24 0 3 6.5  0 1 16.5
14 24 0 2 2.5  0 1 11  
15 25 0 2 3    0 1 12  
16 1  0 0 0    0 1 0   
];

end

Here's how you call the solver:

import GasModels
import HiGHS
fname = "model30ss.m"
slpopt = GasModels._SLP.Optimizer(HiGHS.Optimizer; tol = 1e-5)
result = GasModels.solve_ogf(fname, GasModels.WPGasModel, slpopt)
display(result["solution"]["receipt"])

Here's the output I get:

iter    objective   inf_pr    inf_du   inf_comp
   0    0.000e+00  2.16e+03  1.00e+00  0.00e+00
   1    1.852e+03  8.80e+00  1.19e+01  0.00e+00
   2    1.782e+03  3.05e-01  1.00e+00  5.50e-02
   3    1.782e+03  9.98e-02  1.00e+00  5.50e-02
   4    1.782e+03  1.58e-01  1.00e+00  5.50e-02
   5    1.782e+03  5.50e-02  1.00e+00  5.50e-02
   6    6.400e+02  1.37e+01  0.00e+00  0.00e+00
   7    6.400e+02  3.93e+01  0.00e+00  0.00e+00
   8    6.400e+02  1.42e-13  0.00e+00  0.00e+00
Dict{String, Dict{String, Float64}} with 16 entries:
  "4" => Dict("fg"=>-0.0)
  "1" => Dict("fg"=>-0.0)
  "12" => Dict("fg"=>-0.0)
  "2" => Dict("fg"=>0.24)
  "6" => Dict("fg"=>-0.0)
  "11" => Dict("fg"=>-0.0)
  "13" => Dict("fg"=>-0.0)
  "5" => Dict("fg"=>-0.0)
  "15" => Dict("fg"=>-0.0)
  "16" => Dict("fg"=>1.1)
  "14" => Dict("fg"=>-0.0)
  "7" => Dict("fg"=>-0.0)
  "8" => Dict("fg"=>-0.0)
  "10" => Dict("fg"=>-0.0)
  "9" => Dict("fg"=>0.3)
  "3" => Dict("fg"=>-0.0)

[kaarthiksundar]

Things to do, 1. Trust region update 2. A restoration like algorithm to get to a good starting point if the starting point is infeasible so that the algorithm is robust to bad x0.

[skylerreid]

Just a quick comment: the user-facing functions in GMs are being migrated from run_* to solve_* to maintain similarity across the InfrastructureModels packages. It would probably be a good idea to include that going forward in your SLP solver as well.

[Robbybp]

TODO: I should use InfrastructureModels' solution constructor here

[Robbybp]

This is tricky because IM.build_solution depends on JuMP.value, which uses the MOI.VariablePrimal attribute. I think I need to give the model a mock backend then set this attribute with something like

MOI.set(JuMP.backend(model), MOI.VariablePrimal(1), var, slp.primal_solution[var])

[Robbybp]

I couldn't get this to work in a quick attempt. My best attempt has been:

JuMP.set_optimizer(gm.model, MOIU.MockOptimizer)
MOIU.attach_optimizer(gm.model)
MOI.set(JuMP.unsafe_backend(gm.model), MOI.TerminationStatus(), slp_results.termination_status)
for var in JuMP.all_variables(gm.model)
    val = slp_results.primal_solution[var]
    MOI.set(JuMP.unsafe_backend(gm.model), MOI.VariablePrimal(1), JuMP.index(var), val)
end
gm.model.is_model_dirty = false
solution = _IM.build_solution(
    gm;
    post_processors = [
        sol_psqr_to_p!,
        sol_compressor_p_to_r!,
        sol_regulator_p_to_r!,
    ],
)

This gives me the error message LoadError: No mock primal is set for variable MOI.VariableIndex(12345681), which I don't know how to interpret.

Options I see for solution building:

1. Build solution manually (current approach)
2. Patch InfrastructureModels to accept something like a var_value function or dict to override the use of JuMP.value
3. Eventually make this a proper MOI solver, then update the solve_ogf method to use build_solution

I currently like option 3.