Problem 6.3#
Integrated Energy Grids
Problem 6.3
Assume the gas transport system that is shown in Figure 2. Gas can be injected in Node 1 at a marginal price of 10 EUR/kg or at Node 5 at 15 EUR/kg. There is a demand of 2·\(10^8\) kg/s in Node 3. Finally, Node 2-3 represents a pumping station where gas pressure can be increased by \(k_{23}\)=1.2 at a cost that is proportional to the mass flow \(m_{23}\) by a constant \(o_{23}\) =5 EUR/kg.
To calculate the speed of sound in gas \(c\), assume the universal gas constant \(R\)=8.314 J/molK, the molar mass of methane \(M\)=16 g/mol, compression factor \(Z\)=1.31, and temperature T=25\(^{\circ}\)C. Assume the diameter of the pipelines is \(D\) = 1000 mm, the length \(L_{12}\)=\(L_{34}\)=\(L_{45}\)=100 km, and calculate the Darcy friction coefficient \(f_{D}\) assuming fully-turbulent flow and roughness \(\epsilon\)=0.001m.
(a) Write the equations of the optimization problem that minimizes the total system cost and use the Weymouth equation to represent pressure drop in the pipelines as a function of the mass flow \(m_{ij}\).
We define the pressure at differet node as \(\pi_i\) and the mass flow at different pipes as \(m_{ij}\) We have 9 variables (\(m_{12}\), \(m_{23}\), \(m_{34}\), \(m_{45}\), \(\pi_1\), \(\pi_2\), \(\pi_3\), \(\pi_4\), \(\pi_5\)).
We impose the energy balance constraint in nodes 1, 4 and 5, and the pipeline mass conservation and Weymouth equation in pipelines 12, 34 and 45.
where the coefficients \(a_{ij}=\frac{Lf^Dc^2}{DA^2}\) depend on the physical properties of the pipelines and the methane gas.
(b) Use gurobipy to formulate and solve the problem
We start by calculating the cofficients \(a_{ij}\) from the Weymouth equation
import numpy as np
The cross-sectional area \(A\), and the speed of sound in gas \(c\), can be calculated as
D = 1 # m
Z = 1.31
R = 8.314 # J/molk
M = 0.016 # Kg/mol
T = 273+25 # K
A = np.pi*(D/2)**2
c=np.sqrt(Z*R*T/M)
c
np.float64(450.3900615022494)
The Darcy friction coefficient can be calculated as \(\frac{1}{\sqrt f_D}=-2log_{10}(\frac{\epsilon}{3.7D})\)
epsilon = 0.001 # m
f_D=(1/(-2*np.log10(epsilon/(3.7*D))))**2
f_D
np.float64(0.0196354659355267)
The coefficients from the Weymouth equation can be calculated as \(a_{ij}=\frac{Lf_Dc^2}{DA^2}\)
L = 100*1000 # km -> m
L_12=L
a_12=L_12*f_D*c**2/(D*A**2)
a_34=a_12
a_45=a_12
a_12
np.float64(645712279.9219108)
We import gurobipy and write the problem
import gurobipy as gp
from gurobipy import GRB
k_23 = 1.2
d_4 = 2*100000000 #kg/s
o_1 = 10 #EUR/kg
o_5 = 15 #EUR/kg
o_m23 = 5 #EUR/kg
model = gp.Model("My_quadratic_problem")
g_1 = model.addVar(vtype=GRB.CONTINUOUS, lb=0, name="g_1")
g_5 = model.addVar(vtype=GRB.CONTINUOUS, lb=0, name="g_5")
m_12 = model.addVar(vtype=GRB.CONTINUOUS, lb=0, name="m_12")
m_23 = model.addVar(vtype=GRB.CONTINUOUS, lb=0, name="m_23")
m_34 = model.addVar(vtype=GRB.CONTINUOUS, lb=0, name="m_34")
m_45 = model.addVar(vtype=GRB.CONTINUOUS, lb=-200000000, name="m_45")
abs_m_12 = model.addVar(vtype=GRB.CONTINUOUS, lb=0, name="abs_m_12")
abs_m_34 = model.addVar(vtype=GRB.CONTINUOUS, lb=0, name="abs_m_34")
abs_m_45 = model.addVar(vtype=GRB.CONTINUOUS, lb=0, name="abs_m_45")
p_1 = model.addVar(vtype=GRB.CONTINUOUS, lb=0, name="p_1")
p_2 = model.addVar(vtype=GRB.CONTINUOUS, lb=0, name="p_2")
p_3 = model.addVar(vtype=GRB.CONTINUOUS, lb=0, name="p_3")
p_4 = model.addVar(vtype=GRB.CONTINUOUS, lb=0, name="p_4")
p_5 = model.addVar(vtype=GRB.CONTINUOUS, lb=0, name="p_5")
model.setObjective(o_1*g_1 + o_5*g_5 + o_m23*m_23, GRB.MINIMIZE)
constraint_a = model.addLConstr(g_1 == m_12)
constraint_b = model.addLConstr(-g_5 == m_45)
constraint_c = model.addLConstr(-d_4 == -m_34+m_45)
constraint_d = model.addLConstr(m_12 == m_23)
constraint_e = model.addLConstr(m_12 == m_34)
constraint_f = model.addLConstr(p_3 == k_23*p_2)
# We need to define additional constraints to use absolute values
from gurobipy import abs_
model.addConstr(abs_m_12 == abs_(m_12), name='abs_m_12')
model.addConstr(abs_m_34 == abs_(m_34), name='abs_m_34')
model.addConstr(abs_m_45 == abs_(m_45), name='abs_m_45')
constraint_g = model.addQConstr(a_12*m_12*abs_m_12 == p_1**2-p_2**2)
constraint_h = model.addQConstr(a_34*m_34*abs_m_34 == p_3**2-p_4**2)
constraint_i = model.addQConstr(a_45*m_45*abs_m_45 == p_4**2-p_5**2)
model.optimize()
Restricted license - for non-production use only - expires 2026-11-23
Gurobi Optimizer version 12.0.2 build v12.0.2rc0 (linux64 - "Ubuntu 24.04.2 LTS")
CPU model: AMD EPYC 7763 64-Core Processor, instruction set [SSE2|AVX|AVX2]
Thread count: 2 physical cores, 4 logical processors, using up to 4 threads
Optimize a model with 6 rows, 14 columns and 12 nonzeros
Model fingerprint: 0x2eae29bb
Model has 3 quadratic constraints
Model has 3 simple general constraints
3 ABS
Variable types: 14 continuous, 0 integer (0 binary)
Coefficient statistics:
Matrix range [1e+00, 1e+00]
QMatrix range [1e+00, 6e+08]
Objective range [5e+00, 2e+01]
Bounds range [2e+08, 2e+08]
RHS range [2e+08, 2e+08]
Warning: Quadratic constraints contain large coefficients
Consider reformulating model or setting NumericFocus parameter
to avoid numerical issues.
Presolve removed 5 rows and 8 columns
Presolve time: 0.00s
Presolved: 22 rows, 13 columns, 41 nonzeros
Presolved model has 6 bilinear constraint(s)
Warning: Model contains variables with very large bounds participating
in product terms.
Presolve was not able to compute smaller bounds for these variables.
Consider bounding these variables or reformulating the model.
Solving non-convex MIQCP
Variable types: 13 continuous, 0 integer (0 binary)
Root relaxation: objective 3.000000e+09, 1 iterations, 0.00 seconds (0.00 work units)
Nodes | Current Node | Objective Bounds | Work
Expl Unexpl | Obj Depth IntInf | Incumbent BestBd Gap | It/Node Time
0 0 3.0000e+09 0 1 - 3.0000e+09 - - 0s
0 0 3.0000e+09 0 1 - 3.0000e+09 - - 0s
0 2 3.0000e+09 0 1 - 3.0000e+09 - - 0s
370427 6222 infeasible 633 - 3.0000e+09 - 0.3 5s
H374589 202 3.000000e+09 3.0000e+09 0.00% 0.3 5s
Explored 375636 nodes (128925 simplex iterations) in 5.12 seconds (0.86 work units)
Thread count was 4 (of 4 available processors)
Solution count 1: 3e+09
Optimal solution found (tolerance 1.00e-04)
Best objective 3.000000000000e+09, best bound 3.000000000000e+09, gap 0.0000%
We print the values of the variables in the solution
print(str(g_1.VarName) + ' ' + str(g_1.x))
print(str(g_5.VarName) + ' ' + str(g_5.x))
print(str(m_12.VarName) + ' ' + str(m_12.x))
print(str(m_23.VarName) + ' ' + str(m_23.x))
print(str(m_34.VarName) + ' ' + str(m_34.x))
print(str(m_45.VarName) + ' ' + str(m_45.x))
print(str(p_1.VarName) + ' ' + str(round(p_1.x)) + ' Pa')
print(str(p_2.VarName) + ' ' + str(round(p_2.x)) + ' Pa')
print(str(p_3.VarName) + ' ' + str(round(p_3.x)) + ' Pa')
print(str(p_4.VarName) + ' ' + str(round(p_4.x)) + ' Pa')
print(str(p_5.VarName) + ' ' + str(round(p_5.x)) + ' Pa')
g_1 0.0021958809035774597
g_5 199999999.99780396
m_12 0.0021958809035774597
m_23 0.0021958809035774597
m_34 0.0021958809035774597
m_45 -199999999.99780396
p_1 73 Pa
p_2 46 Pa
p_3 56 Pa
p_4 0 Pa
p_5 5082173865218 Pa