Scippy

SCIP

Solving Constraint Integer Programs

gastrans.c
Go to the documentation of this file.
1 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
2 /* */
3 /* This file is part of the program and library */
4 /* SCIP --- Solving Constraint Integer Programs */
5 /* */
6 /* Copyright 2002-2022 Zuse Institute Berlin */
7 /* */
8 /* Licensed under the Apache License, Version 2.0 (the "License"); */
9 /* you may not use this file except in compliance with the License. */
10 /* You may obtain a copy of the License at */
11 /* */
12 /* http://www.apache.org/licenses/LICENSE-2.0 */
13 /* */
14 /* Unless required by applicable law or agreed to in writing, software */
15 /* distributed under the License is distributed on an "AS IS" BASIS, */
16 /* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. */
17 /* See the License for the specific language governing permissions and */
18 /* limitations under the License. */
19 /* */
20 /* You should have received a copy of the Apache-2.0 license */
21 /* along with SCIP; see the file LICENSE. If not visit scipopt.org. */
22 /* */
23 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
24 
25 /**@file gastrans.c
26  * @brief Simple Gas Transportation Model
27  * @author Stefan Vigerske
28  *
29  * This example shows how to setup nonlinear constraints with signpower-expression when using SCIP as callable library.
30  * The example implements a model for the distribution of gas through a network of pipelines, which
31  * is formulated as a cost minimization subject to nonlinear flow-pressure relations, material balances,
32  * and pressure bounds. The Belgian gas network is used as an example.
33  *
34  * The model is taken from the GAMS model library:
35  * https://www.gams.com/latest/gamslib_ml/libhtml/gamslib_gastrans.html
36  *
37  * Original model source:
38  * @par
39  * D. de Wolf and Y. Smeers@n
40  * The Gas Transmission Problem Solved by and Extension of the Simplex Algorithm@n
41  * Management Science 46, 11 (2000), 1454-1465
42  */
43 
44 /*--+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8----+----9----+----0----+----1----+----2*/
45 
46 #include <stdio.h>
47 #include <math.h>
48 
49 #include "scip/scip.h"
50 #include "scip/scipdefplugins.h"
51 
52 /* Model parameters */
53 
54 /** node data structure */
55 typedef struct NodeData {
56  const char* name;
62 } NodeData;
63 
64 /** arc data structure */
65 typedef struct ArcData {
66  int node1;
67  int node2;
71 } ArcData;
72 
73 /** number of nodes (towns) */
74 #define nnodes 20
75 
76 /** number of arcs */
77 #define narcs 24
78 
79 /** value we use to represent infinity */
80 #define infinity 1e+20
81 
82 /** data of nodes */
83 static const NodeData nodedata[] =
84 {
85  /* name supplylo supplyup pressurelo pressureup cost */
86  {"Anderlues", 0.0, 1.2, 0.0, 66.2, 0.0 }, /* 0 */
87  {"Antwerpen", -infinity, -4.034, 30.0, 80.0, 0.0 }, /* 1 */
88  {"Arlon", -infinity, -0.222, 0.0, 66.2, 0.0 }, /* 2 */
89  {"Berneau", 0.0, 0.0, 0.0, 66.2, 0.0 }, /* 3 */
90  {"Blaregnies", -infinity, -15.616, 50.0, 66.2, 0.0 }, /* 4 */
91  {"Brugge", -infinity, -3.918, 30.0, 80.0, 0.0 }, /* 5 */
92  {"Dudzele", 0.0, 8.4, 0.0, 77.0, 2.28 }, /* 6 */
93  {"Gent", -infinity, -5.256, 30.0, 80.0, 0.0 }, /* 7 */
94  {"Liege", -infinity, -6.385, 30.0, 66.2, 0.0 }, /* 8 */
95  {"Loenhout", 0.0, 4.8, 0.0, 77.0, 2.28 }, /* 9 */
96  {"Mons", -infinity, -6.848, 0.0, 66.2, 0.0 }, /* 10 */
97  {"Namur", -infinity, -2.120, 0.0, 66.2, 0.0 }, /* 11 */
98  {"Petange", -infinity, -1.919, 25.0, 66.2, 0.0 }, /* 12 */
99  {"Peronnes", 0.0, 0.96, 0.0, 66.2, 1.68 }, /* 13 */
100  {"Sinsin", 0.0, 0.0, 0.0, 63.0, 0.0 }, /* 14 */
101  {"Voeren", 20.344, 22.012, 50.0, 66.2, 1.68 }, /* 15 */
102  {"Wanze", 0.0, 0.0, 0.0, 66.2, 0.0 }, /* 16 */
103  {"Warnand", 0.0, 0.0, 0.0, 66.2, 0.0 }, /* 17 */
104  {"Zeebrugge", 8.87, 11.594, 0.0, 77.0, 2.28 }, /* 18 */
105  {"Zomergem", 0.0, 0.0, 0.0, 80.0, 0.0 } /* 19 */
106 };
107 
108 /** data of arcs */
109 static const ArcData arcdata[] =
110 {
111  /* node1 node2 diameter length active */
112  { 18, 6, 890.0, 4.0, FALSE },
113  { 18, 6, 890.0, 4.0, FALSE },
114  { 6, 5, 890.0, 6.0, FALSE },
115  { 6, 5, 890.0, 6.0, FALSE },
116  { 5, 19, 890.0, 26.0, FALSE },
117  { 9, 1, 590.1, 43.0, FALSE },
118  { 1, 7, 590.1, 29.0, FALSE },
119  { 7, 19, 590.1, 19.0, FALSE },
120  { 19, 13, 890.0, 55.0, FALSE },
121  { 15, 3, 890.0, 5.0, TRUE },
122  { 15, 3, 395.0, 5.0, TRUE },
123  { 3, 8, 890.0, 20.0, FALSE },
124  { 3, 8, 395.0, 20.0, FALSE },
125  { 8, 17, 890.0, 25.0, FALSE },
126  { 8, 17, 395.0, 25.0, FALSE },
127  { 17, 11, 890.0, 42.0, FALSE },
128  { 11, 0, 890.0, 40.0, FALSE },
129  { 0, 13, 890.0, 5.0, FALSE },
130  { 13, 10, 890.0, 10.0, FALSE },
131  { 10, 4, 890.0, 25.0, FALSE },
132  { 17, 16, 395.5, 10.5, FALSE },
133  { 16, 14, 315.5, 26.0, TRUE },
134  { 14, 2, 315.5, 98.0, FALSE },
135  { 2, 12, 315.5, 6.0, FALSE }
136 };
137 
138 /** gas temperatur (K) */
139 static const SCIP_Real gastemp = 281.15;
140 
141 /** absolute rugosity (mm) */
142 static const SCIP_Real rugosity = 0.05;
143 
144 /** density of gas relative to air */
145 static const SCIP_Real density = 0.616;
146 
147 /** compressibility factor */
148 static const SCIP_Real compressibility = 0.8;
149 
150 
151 
152 /** sets up problem */
153 static
155  SCIP* scip /**< SCIP data structure */
156  )
157 {
158  SCIP_VAR* flow[narcs];
159  SCIP_VAR* supply[nnodes];
160  SCIP_VAR* pressure[nnodes];
161  SCIP_VAR* pressurediff[narcs];
162 
163  SCIP_CONS* flowbalance[nnodes];
164  SCIP_CONS* pressurediffcons[narcs];
165  SCIP_CONS* pressureloss[narcs];
166 
167  char name[SCIP_MAXSTRLEN];
168  int i;
169  int j;
170 
171  /* create empty problem */
172  SCIP_CALL( SCIPcreateProbBasic(scip, "gastrans") );
173 
174  /* create variables for flows and add to problem */
175  for( i = 0; i < narcs; ++i )
176  {
177  (void) SCIPsnprintf(name, SCIP_MAXSTRLEN, "flow_%s_%s", nodedata[arcdata[i].node1].name, nodedata[arcdata[i].node2].name);
178  SCIP_CALL( SCIPcreateVarBasic(scip, &flow[i], name, arcdata[i].active ? 0.0 : -SCIPinfinity(scip), SCIPinfinity(scip), 0.0, SCIP_VARTYPE_CONTINUOUS) );
179 
180  SCIP_CALL( SCIPaddVar(scip, flow[i]) );
181  }
182 
183  /* create variables for pressure difference and add to problem */
184  for( i = 0; i < narcs; ++i )
185  {
186  (void) SCIPsnprintf(name, SCIP_MAXSTRLEN, "pressurediff_%s_%s", nodedata[arcdata[i].node1].name, nodedata[arcdata[i].node2].name);
187  SCIP_CALL( SCIPcreateVarBasic(scip, &pressurediff[i], name, -SCIPinfinity(scip), SCIPinfinity(scip), 0.0, SCIP_VARTYPE_CONTINUOUS) );
188 
189  SCIP_CALL( SCIPaddVar(scip, pressurediff[i]) );
190  }
191 
192  /* create variables for supply at nodes and add to problem */
193  for( i = 0; i < nnodes; ++i )
194  {
195  (void) SCIPsnprintf(name, SCIP_MAXSTRLEN, "supply_%s", nodedata[i].name);
196  SCIP_CALL( SCIPcreateVarBasic(scip, &supply[i], name,
197  nodedata[i].supplylower == -infinity ? -SCIPinfinity(scip) : nodedata[i].supplylower,
198  nodedata[i].supplyupper == infinity ? SCIPinfinity(scip) : nodedata[i].supplyupper,
199  nodedata[i].cost, SCIP_VARTYPE_CONTINUOUS) ); /*lint !e777*/
200 
201  SCIP_CALL( SCIPaddVar(scip, supply[i]) );
202  }
203 
204  /* create variables for squared pressure at nodes and add to problem */
205  for( i = 0; i < nnodes; ++i )
206  {
207  (void) SCIPsnprintf(name, SCIP_MAXSTRLEN, "pressure_%s", nodedata[i].name);
208  SCIP_CALL( SCIPcreateVarBasic(scip, &pressure[i], name,
209  nodedata[i].pressurelower*nodedata[i].pressurelower, nodedata[i].pressureupper*nodedata[i].pressureupper,
210  0.0, SCIP_VARTYPE_CONTINUOUS) );
211 
212  SCIP_CALL( SCIPaddVar(scip, pressure[i]) );
213  }
214 
215  /* create flow balance constraints and add to problem
216  * for each node i: outflows - inflows = supply
217  */
218  for( i = 0; i < nnodes; ++i )
219  {
220  SCIP_Real minusone;
221 
222  minusone = -1.0;
223 
224  (void) SCIPsnprintf(name, SCIP_MAXSTRLEN, "flowbalance%s", nodedata[i].name);
225  SCIP_CALL( SCIPcreateConsBasicLinear(scip, &flowbalance[i], name, 1, &supply[i], &minusone, 0.0, 0.0) );
226 
227  for( j = 0; j < narcs; ++j )
228  {
229  /* check for outflows */
230  if( arcdata[j].node1 == i )
231  {
232  SCIP_CALL( SCIPaddCoefLinear(scip, flowbalance[i], flow[j], +1.0) );
233  }
234 
235  /* check for inflows */
236  if( arcdata[j].node2 == i )
237  {
238  SCIP_CALL( SCIPaddCoefLinear(scip, flowbalance[i], flow[j], -1.0) );
239  }
240  }
241 
242  SCIP_CALL( SCIPaddCons(scip, flowbalance[i]) );
243  }
244 
245  /* create pressure difference constraints and add to problem
246  * pressurediff[node1 to node2] = pressure[node1] - pressure[node2]
247  */
248  for( i = 0; i < narcs; ++i )
249  {
250  SCIP_VAR* vars[3];
251  SCIP_Real coefs[3] = { 1.0, -1.0, 1.0 };
252 
253  assert(pressurediff[i] != NULL);
254  assert(arcdata[i].node1 >= 0 && arcdata[i].node1 < nnodes);
255  assert(arcdata[i].node2 >= 0 && arcdata[i].node2 < nnodes);
256  assert(pressure[arcdata[i].node1] != NULL);
257  assert(pressure[arcdata[i].node2] != NULL);
258 
259  vars[0] = pressurediff[i];
260  vars[1] = pressure[arcdata[i].node1];
261  vars[2] = pressure[arcdata[i].node2];
262 
263  (void) SCIPsnprintf(name, SCIP_MAXSTRLEN, "pressurediff_%s_%s", nodedata[arcdata[i].node1].name, nodedata[arcdata[i].node2].name);
264  SCIP_CALL( SCIPcreateConsBasicLinear(scip, &pressurediffcons[i], name, 3, vars, coefs, 0.0, 0.0) );
265 
266  SCIP_CALL( SCIPaddCons(scip, pressurediffcons[i]) );
267  }
268 
269  /* create pressure loss constraints and add to problem
270  * pressure loss on active arcs: flow(i) * flow(i) + coef * pressurediff(i) <= 0.0
271  * on regular pipes: flow(i) * abs(flow(i)) - coef * pressurediff(i) == 0.0,
272  * where coef = 96.074830e-15*power(diameter(i)^5/lambda/compressibility/temperatur/length(i)/density
273  * and lambda = (2*log10(3.7*diameter(i)/rugosity))^(-2);
274  */
275  for( i = 0; i < narcs; ++i )
276  {
277  SCIP_EXPR* exprflow;
278  SCIP_EXPR* exprs[2];
279  SCIP_Real coefs[2];
280  SCIP_EXPR* exprsum;
281  SCIP_Real coef;
282 
283  coef = 96.074830e-15 * pow(arcdata[i].diameter, 5.0) * pow(2.0*log10(3.7*arcdata[i].diameter / rugosity), 2.0)
284  / compressibility / gastemp / arcdata[i].length / density;
285 
286  /* we can always use the signpower-expression, because flow(i) is positive for active arcs */
287  SCIP_CALL( SCIPcreateExprVar(scip, &exprflow, flow[i], NULL, NULL) );
288  SCIP_CALL( SCIPcreateExprSignpower(scip, &exprs[0], exprflow, 2.0, NULL, NULL) );
289  SCIP_CALL( SCIPreleaseExpr(scip, &exprflow) );
290 
291  SCIP_CALL( SCIPcreateExprVar(scip, &exprs[1], pressurediff[i], NULL, NULL) );
292 
293  coefs[0] = 1.0;
294  coefs[1] = arcdata[i].active ? coef : -coef;
295  SCIP_CALL( SCIPcreateExprSum(scip, &exprsum, 2, exprs, coefs, 0.0, NULL, NULL) );
296  SCIP_CALL( SCIPreleaseExpr(scip, &exprs[1]) );
297  SCIP_CALL( SCIPreleaseExpr(scip, &exprs[0]) );
298 
299  (void) SCIPsnprintf(name, SCIP_MAXSTRLEN, "pressureloss_%s_%s", nodedata[arcdata[i].node1].name, nodedata[arcdata[i].node2].name);
300  SCIP_CALL( SCIPcreateConsBasicNonlinear(scip, &pressureloss[i], name, exprsum, arcdata[i].active ? -SCIPinfinity(scip) : 0.0, 0.0) );
301  SCIP_CALL( SCIPreleaseExpr(scip, &exprsum) );
302 
303  SCIP_CALL( SCIPaddCons(scip, pressureloss[i]) );
304  }
305 
306  /* release variables and constraints
307  * the problem has them captured, and we do not require them anymore
308  */
309  for( i = 0; i < nnodes; ++i )
310  {
311  SCIP_CALL( SCIPreleaseVar(scip, &supply[i]) );
312  SCIP_CALL( SCIPreleaseVar(scip, &pressure[i]) );
313 
314  SCIP_CALL( SCIPreleaseCons(scip, &flowbalance[i]) );
315  }
316 
317  for( i = 0; i < narcs; ++i )
318  {
319  SCIP_CALL( SCIPreleaseVar(scip, &flow[i]) );
320  SCIP_CALL( SCIPreleaseVar(scip, &pressurediff[i]) );
321 
322  SCIP_CALL( SCIPreleaseCons(scip, &pressurediffcons[i]) );
323  SCIP_CALL( SCIPreleaseCons(scip, &pressureloss[i]) );
324  }
325 
326  return SCIP_OKAY;
327 }
328 
329 /** runs gas transportation example */
330 static
332 {
333  SCIP* scip;
334 
335  SCIP_CALL( SCIPcreate(&scip) );
337 
338  SCIPinfoMessage(scip, NULL, "\n");
339  SCIPinfoMessage(scip, NULL, "*****************************************\n");
340  SCIPinfoMessage(scip, NULL, "* Running Gas Transportation Model *\n");
341  SCIPinfoMessage(scip, NULL, "*****************************************\n");
342  SCIPinfoMessage(scip, NULL, "\n");
343 
344  SCIP_CALL( setupProblem(scip) );
345 
346  SCIPinfoMessage(scip, NULL, "Original problem:\n");
347  SCIP_CALL( SCIPprintOrigProblem(scip, NULL, "cip", FALSE) );
348 
349  SCIPinfoMessage(scip, NULL, "\n");
350  SCIP_CALL( SCIPpresolve(scip) );
351 
352  /* SCIPinfoMessage(scip, NULL, "Presolved problem:\n");
353  SCIP_CALL( SCIPprintTransProblem(scip, NULL, "cip", FALSE) );
354  */
355 
356  SCIPinfoMessage(scip, NULL, "\nSolving...\n");
357  SCIP_CALL( SCIPsolve(scip) );
358 
359  /*
360  if( SCIPgetNSols(scip) > 0 )
361  {
362  SCIPinfoMessage(scip, NULL, "\nSolution:\n");
363  SCIP_CALL( SCIPprintSol(scip, SCIPgetBestSol(scip), NULL, FALSE) );
364  }
365  */
366 
367  SCIP_CALL( SCIPfree(&scip) );
368 
369  return SCIP_OKAY;
370 }
371 
372 
373 /** main method starting SCIP */
374 int main(
375  int argc, /**< number of arguments from the shell */
376  char** argv /**< array of shell arguments */
377  )
378 { /*lint --e{715}*/
379  SCIP_RETCODE retcode;
380 
381  retcode = runGastrans();
382 
383  /* evaluate return code of the SCIP process */
384  if( retcode != SCIP_OKAY )
385  {
386  /* write error back trace */
387  SCIPprintError(retcode);
388  return -1;
389  }
390 
391  return 0;
392 }
#define narcs
Definition: gastrans.c:77
SCIP_Real cost
Definition: gastrans.c:61
SCIP_Real supplyupper
Definition: gastrans.c:58
SCIP_RETCODE SCIPcreateConsBasicLinear(SCIP *scip, SCIP_CONS **cons, const char *name, int nvars, SCIP_VAR **vars, SCIP_Real *vals, SCIP_Real lhs, SCIP_Real rhs)
#define infinity
Definition: gastrans.c:80
SCIP_Real pressureupper
Definition: gastrans.c:60
#define SCIP_MAXSTRLEN
Definition: def.h:302
SCIP_RETCODE SCIPcreateExprSignpower(SCIP *scip, SCIP_EXPR **expr, SCIP_EXPR *child, SCIP_Real exponent, SCIP_DECL_EXPR_OWNERCREATE((*ownercreate)), void *ownercreatedata)
Definition: expr_pow.c:3199
SCIP_Real supplylower
Definition: gastrans.c:57
SCIP_RETCODE SCIPreleaseVar(SCIP *scip, SCIP_VAR **var)
Definition: scip_var.c:1254
static SCIP_RETCODE runGastrans(void)
Definition: gastrans.c:331
#define FALSE
Definition: def.h:96
SCIP_Real pressurelower
Definition: gastrans.c:59
SCIP_Real length
Definition: gastrans.c:69
SCIP_Real SCIPinfinity(SCIP *scip)
int SCIPsnprintf(char *t, int len, const char *s,...)
Definition: misc.c:10764
#define TRUE
Definition: def.h:95
static const SCIP_Real gastemp
Definition: gastrans.c:139
enum SCIP_Retcode SCIP_RETCODE
Definition: type_retcode.h:63
SCIP_RETCODE SCIPcreateVarBasic(SCIP *scip, SCIP_VAR **var, const char *name, SCIP_Real lb, SCIP_Real ub, SCIP_Real obj, SCIP_VARTYPE vartype)
Definition: scip_var.c:194
static GRAPHNODE ** active
struct NodeData NodeData
SCIP_RETCODE SCIPcreateExprVar(SCIP *scip, SCIP_EXPR **expr, SCIP_VAR *var, SCIP_DECL_EXPR_OWNERCREATE((*ownercreate)), void *ownercreatedata)
Definition: expr_var.c:390
SCIP_RETCODE SCIPcreate(SCIP **scip)
Definition: scip_general.c:292
SCIP_RETCODE SCIPcreateExprSum(SCIP *scip, SCIP_EXPR **expr, int nchildren, SCIP_EXPR **children, SCIP_Real *coefficients, SCIP_Real constant, SCIP_DECL_EXPR_OWNERCREATE((*ownercreate)), void *ownercreatedata)
Definition: expr_sum.c:1079
SCIP_RETCODE SCIPaddCoefLinear(SCIP *scip, SCIP_CONS *cons, SCIP_VAR *var, SCIP_Real val)
void SCIPinfoMessage(SCIP *scip, FILE *file, const char *formatstr,...)
Definition: scip_message.c:208
SCIP_RETCODE SCIPcreateProbBasic(SCIP *scip, const char *name)
Definition: scip_prob.c:180
static const NodeData nodedata[]
Definition: gastrans.c:83
SCIP_RETCODE SCIPsolve(SCIP *scip)
Definition: scip_solve.c:2622
SCIP_RETCODE SCIPaddCons(SCIP *scip, SCIP_CONS *cons)
Definition: scip_prob.c:2778
SCIP_Bool active
Definition: gastrans.c:70
static const SCIP_Real rugosity
Definition: gastrans.c:142
SCIP_RETCODE SCIPcreateConsBasicNonlinear(SCIP *scip, SCIP_CONS **cons, const char *name, SCIP_EXPR *expr, SCIP_Real lhs, SCIP_Real rhs)
SCIP_RETCODE SCIPpresolve(SCIP *scip)
Definition: scip_solve.c:2452
SCIP_RETCODE SCIPprintOrigProblem(SCIP *scip, FILE *file, const char *extension, SCIP_Bool genericnames)
#define NULL
Definition: lpi_spx1.cpp:164
#define SCIP_CALL(x)
Definition: def.h:393
SCIP_Real diameter
Definition: gastrans.c:68
int node2
Definition: gastrans.c:67
#define SCIP_Bool
Definition: def.h:93
SCIP_RETCODE SCIPincludeDefaultPlugins(SCIP *scip)
SCIP_RETCODE SCIPreleaseExpr(SCIP *scip, SCIP_EXPR **expr)
Definition: scip_expr.c:1416
struct ArcData ArcData
static const ArcData arcdata[]
Definition: gastrans.c:109
static SCIP_RETCODE setupProblem(SCIP *scip)
Definition: gastrans.c:154
const char * name
Definition: gastrans.c:56
static const SCIP_Real density
Definition: gastrans.c:145
int node1
Definition: gastrans.c:66
int main(int argc, char **argv)
Definition: gastrans.c:374
SCIP_RETCODE SCIPaddVar(SCIP *scip, SCIP_VAR *var)
Definition: scip_prob.c:1676
SCIP_RETCODE SCIPreleaseCons(SCIP *scip, SCIP_CONS **cons)
Definition: scip_cons.c:1119
#define SCIP_Real
Definition: def.h:186
#define nnodes
Definition: gastrans.c:74
void SCIPprintError(SCIP_RETCODE retcode)
Definition: scip_general.c:220
default SCIP plugins
SCIP callable library.
SCIP_RETCODE SCIPfree(SCIP **scip)
Definition: scip_general.c:324
static const SCIP_Real compressibility
Definition: gastrans.c:148