Optimal control of compressor stations in a coupled gas-to-power network

We introduce a tool for simulation and optimization of gas pipeline networks coupled to power grids by gas-to-power plants. The model under consideration consists of the isentropic Euler equations to describe the gas flow coupled to the AC powerflow equations. A compressor station is installed to control the gas pressure such that certain bounds are satisfied. A numerical case study is presented that showcases effects of fast changes in power demand on gas pipelines and necessary operator actions.

because these plants are relatively flexible in comparison to coal or nuclear plants. A welcome advantage of this approach is the possibility to run gas plants with fuel produced from renewable electricity via power-togas plants, thereby reducing or even negating carbon emissions of the gas plants. For a review of power-togas capability see [6].
It is desirable to have a joint optimal control framework for power and gas sector of the energy system to model this compensation. So far only steady-state flow in the gas network has been considered [3,18,20], which may be too coarse for several applications. Therefore, we focus on an optimal control strategy for the instationary gas network model [1] coupled to a power grid [2] via compressor stations, see for example [8,15,17]. The mathematical foundation for the gas-to-power coupling has been recently introduced in [5], where conditions for the well-posedness have been derived and proved. This work is the first attempt to model this interaction and yields an understanding of the underlying equations. The next step will be real-world scenarios.

Optimal control problem
The gas dynamics within each pipeline of the considered gas networks are modeled by the isentropic Euler equations, supplemented with suitable coupling and boundary conditions. For the power grid, we apply the wellknown powerflow equations. The coupling between gas and power networks at gas-driven power plants is modeled by (algebraic) demand-dependent gas consumption terms. To react on the demand-dependent influences on the gas network, controllable devices as compressor stations are considered within the gas network. The aim is to fulfill given state restrictions like pressure bounds whereas at the same time the entire fuel gas or power arXiv:1901.00522v2 [math.OC] 22 Mar 2019 consumption of the compressor stations is to be minimized.
The network under consideration is similar to the one presenetd in [5] and is depicted in Figure 1. In addition, there is now a compressor with a time-dependent control u(t), which is used to satisfy pressure bounds. The optimal control problem is to minimize compressor costs while satisfying power demand and gas dynamics: Mathematically, this is an instationary nonlinear optimization problem constrained by partial differential equations, see [9] for an overview. To solve the problem, we make use of a first-discretize-then-optimize approach and apply the interior point solver IPOPT [16]. The necessary gradient information for IPOPT, i.e., gradients with respect to all controllable devices, is efficiently computed via adjoint equations. Here, the underlying systems can be solved time-step-wise (backwards in time), where additionally the sparsity structure is exploited.
We remark that the cost function is given in Subsection 2.3, while the bounds on the pressure are introduced as box constraints within the numerical optimization procedure in Section 3. Within the Subsections 2.1-2.5, we now describe the constraints of the optimal control problem in detail and focus on the technical details.

The isentropic Euler equations
The gas network is modeled by the isentropic Euler equations (see [4,1]), which govern gas flow in each pipeline between nodes, where ρ is the density, q = ρv is the flow, p(ρ) = κρ γ is the pressure function. In our example we use γ = 1 and κ = c 2 = (340 m s ) 2 , that is, the isothermal Euler equations with speed of sound c. The Euler equations must be met for 0 ≤ x ≤ l e and t ≥ 0, where l e is the length of the pipe. Furthermore, S is a source term, S0 S17 S4 S8 S20 S25 Fig. 1: Gas network connected to a power grid. Red nodes are PQ/demand nodes, green nodes are generators (PV nodes) and the blue node is the slack bus (also a generator, with gas consumption of the form ε(P ) = a 0 + a 1 P + a 2 P 2 ). The circle symbol indicates a compressor station.
where λ is the solution to the Prandtl-Colebrook formula, The Reynolds number Re is given by with dynamic viscosity η = 10 −5 kg ms .
The roughness k e and diameter d e of the pipes are all the same in our example and given by k e = 5 · 10 −4 m, d e = 6 · 10 −1 m.

Coupling at gas nodes
At the nodes we use the usual Kirchhoff-type coupling conditions: The pressure is the same near the node in all pipes connected to it and the flows must add up to zero (where the sign for inflow is positive, the sign for outflow negative), ingoing pipes The example gas network we are using is a small part of the GasLib-40 network [10]. In Table 1 the only remaining parameter, the length l e of each pipe, is gathered. Note that no length is given for the arc connecting S0 and S17, because only a compressor is situated between these nodes.

Compressor stations
To compensate for pressure losses in the gas network, we consider compressor stations, which are also modelled as arcs. Those arcs have (time-dependent) in-and outgoing pressure (p in , p out ) and flux values (q in , q out ). In general, the two separate flux values allow the modelling of fuel gas consumption of the compressor station, whereas we will consider an external power supply for the compressor and therefore have q in = q out . The power consumption is modelled as a quadratic function of the power required for the compression process. Denoting this as function P (p in (t), q in (t), p out (t), q out (t)), our objective function in the optimal control problem is of the form Note that the power consumption does not influence the network dynamics and is therefore only of interest for the optimization procedure. For our investigations below it is sufficient to know that the consumption and therewith the costs increase if the ratio p out /p in increases. For the details of the power consumption model, we refer to [17, Section 3.2.3]. The influence of the compressor station on the network dynamics is modelled by the control u(t) of the pressure difference:

Power model
For the power grid we use the AC powerflow equations (see [7] for an introduction), where P k , Q k are real and reactive power at node k,|V k | is the voltage amplitude, φ k is the phase (and φ k,j = φ k − φ j ) and B kj , G kj are parameters of the transmission lines between nodes k and j or of the node k for B kk and G kk . Each node is either the slack bus (in our case N1; V and φ given), a generator bus (N2 and N3; V and P given) or a load bus (N4 through N9; P and Q given). All in all for N nodes we get 2N equations for 2N variables.
The considered power grid is taken from the example "case9" of the MATPOWER Matlab programming suite [19]. A per-unit system is used, whose base power and voltage are 100MW and 345kV respectively. The corresponding node and transmission line parameters are gathered in Table 2, these are the entries of the nodal admittance matrix (see [7]).

Coupling
The last ingredient is a model for converting gas to power at a gas power plant. In our example, it will be situated between the nodes S4 and N1 and convert a gas flow ε into a real power output P according to ε(P ) = a 0 + a 1 P + a 2 P 2 .
The flow ε must be diverted from the pipeline network, hence the coupling condition at node S4 must be changed to The details of simulation of such a combined network and the treatment of all arising mathematical issues can be found in [5]. We now showcase a concrete example.

Problem setup
As already noted, we consider a small part of the GasLib-40 network from [2] consisting of 7 pipelines with a total length of 152km. This network is extended by a compressor station and additionally connected to a power grid with 9 nodes by a gas-to-power generator. For this coupled gas-power network, we simulate a sudden increase in power demand within the power grid and study its effect on the gas network. The considered compressor station is supposed to compensate part of the pressure losses in the gas network such that a given pressure bound is satisfied all the time, while power consumption of the compressor is minimized.
To complete the problem description, the following initial and boundary conditions are given: -P (t) and Q(t) at load nodes, -P (t) and V (t) at generator nodes, -V (t) and φ(t) at the slack bus, -p(t) at S5, -q(t) at S25, -p(x, 0), q(x, 0) for all pipelines.
More precisely, the initial conditions for the power network are given in Table 3. These remain constant over time except for the power demand at node N5, which changes linearly between t = 1 hour and t = 1.5 hours from 0.9 p.u. to 1.8 p.u. for the real power and from 0.3 p.u. to 0.6 p.u. for the reactive power, see also Figure 2.
For the gas network the incoming pressure at S5 is fixed at 60bar, the outflow at S25 is fixed at q = Table 3: Initial conditions of the power grid (p.u.) The fuel consumption parameters we use in equation (5) are given by a 0 = 2, a 1 = 5, a 2 = 10. Since the data for the considered gas and power networks are taken from different sources, the parameters of the gas-to-power generator are chosen in such a way that a significant influence is caused.
Further, the pressure at S25 is supposed to satisfy a lower pressure bound of 41bar, i.e., for all times t, where we consider a time horizon of T = 12 hours.
As we will see below, the compressor station between nodes S0 and S17 will have to run at a certain time, i.e., u(t) > 0, to keep this pressure bound. In general, a solution to the described optimal control problem consists of the control u(t) and the entire network state for all times t, i.e., -P (t), Q(t), V (t), φ(t) for all nodes in the power grid, -p(x, t), q(x, t) for all pipelines in the gas network, fulfilling the given model equations and the pressure constraint.

Discretization and Optimization Schemes
To solve the described optimal control problem, we follow a first-discretize-then-optimize approach. The model equations of the power grid only require a discretization in time, which means that the given boundary conditions and the powerflow equations (4) hold for discrete times t j = j∆t with ∆t = 15 minutes and j ∈ {0, . . . , 48} in our scenario. The discretization of the isentropic Euler equations within the pipelines of the gas network additionally requires a spatial grid (here with grid sizes ∆x e ≈ 1 km) and an appropriate discretization scheme. Here we apply an implicit box scheme [14], which allows the application of large time steps as ∆t = 15 minutes for the considered spatial grid.
Considering the isentropic Euler equations as a system of balance laws of the form the applied scheme reads The numerical approximation of the balance law is thought in the following sense: Together with the algebraic equations modelling the compressor station and the coupling and boundary conditions, the discretization process results in a system of nonlinear equations for all state variables of the coupled gas-power network. For simulation purposes, the entire discretized system is solved with Newton's method. Note that the system can be solved time-step per timestep and that we exploit the sparsity structure of the underlying Jacobian matrices.
So far we can only compute the state of the considered gas-power network (for a given compressor control u(t)) and evaluate quantities of interest like the power consumption of the compressor station or the pressure constraint within the time horizon of the simulation. For the given time discretization, the compressor costs (3) are approximated by the trapezoidal rule and formally contained in J(u, y(u)) below. Next, we want to solve the (discretized) optimal control problem, i.e., to find control values u(t j ) such that the pressure constraint is satisfied, while the power consumption of the compressor is minimized. For this purpose, we apply the interior point optimization code IPOPT [16] to the (reduced) discretized optimal control problem. In addition to our simulation procedure to evaluate quantities of interest for a given control, IPOPT further requires gradient information of those quantities with respect to the control. Based on the considered discretization, such information can be efficiently computed by an adjoint approach, which we briefly describe in the following. Therefore, we consider the discretized optimal control problem in the following abstract form: min J(u, y(u)) s.t. u ∈ R Nu , y ∈ R Ny , model equations: E(u, y(u)) = 0, constraints: g(u, y) ≥ 0.
The vector u contains all control variables (here the compressor control u(t j ) for all times t j = j∆t ∈ [0, 24]) and y contains all state variables of the coupled gaspower network of all time steps t j . The mapping u → y(u) is implicitly given by our simulation procedure. Thus, IPOPT does not have to care about the model equations formally summarized in E(u, y), but only about the further constraints g(u, y) and minimizing the objective function J(u, y). Accordingly, we need to provide total derivatives of J and g with respect to u.
In the following, we consider the computation of these derivatives only for the objective function, since the procedure is identical for the constraints, and we follow the description given in [12]. First of all, the chain rule yields While the partial derivatives of J with respect to u and y can be directly computed, the derivatives of the states y with respect to the control u are only implicitly given. Differentiating the model equations E(u, y(u)) = 0 yields 0 = d d u E(u, y(u)) = ∂ ∂u E(u, y(u))+ ∂ ∂y E(u, y(u)) d d u y(u) and therewith (formally) Even though the partial derivatives on the right-handside of (8) can be directly computed, one would have to solve Nu systems of linear equations here. Instead of that, we insert (8) into (7) and get With the so-called adjoint state ξ as the solution of the adjoint equation we finally have d d u J(u, y(u)) = ∂ ∂u J(u, y(u)) + ξ T ∂ ∂u E(u, y(u)) . (10) It is the fact that (9) is a single linear system and has a special structure, which can be easily exploited (see for instance [11,13]), which makes the computation of derivatives via the presented adjoint approach very efficient. Nevertheless note that for given control variables u one still has to solve the model equations to get y(u), before one may compute the gradient information via (9) and (10).

Results
We first discuss the simulation with inactive compressor. In the course of the simulation, due to the increase in power demand at node N5, the power demand at the slack bus rises as well and leads to increased fuel demand at node S4. This increases the inflow into the gas network, as can be seen in Figure 3. Also the pressure in the final node S25 decreases and violates the pressure bound after approximately 4 hours, see Figure 4. After the optimization procedure, the compressor station compensates part of the pressure losses in the gas network such that the pressure bound is satisfied all the time. Since the power consumption of the compressor station is minimized within the optimization, the pressure constraint is active after roughly 4 hours (see again Figure 4), i.e., the compressor station applies as little as possible energy.

Conclusion
The proposed optimization model allows to predict pressure transgressions within a coupled gas-to-power framework. Simulation and optimization tasks are efficiently solved by exploiting the underlying nonlinear problem structure while keeping the full transient regime. This makes it possible to track bounds much more accurately than with a steady state model, thereby achieving lower costs.