Optimal Valve Closure Operations for Pressure Suppression in Fluid Transport Pipelines

I. Introduction

HYDRAULIC shock phenomenon is a pressure surge or shock which commonly occurs when pipeline valves are manipulated in large-scale facilities for fluid distributions [1]. It generates pressure waves that propagate in the pipeline system at a high speed with a significant pressure oscillation, inducing harmful effects on the hydraulic facilities (examples of pipeline damage are shown in Fig. 1). Numerous pipeline collapses have occurred, such as pipelines leakage [2] and equipment failure of spacecraft propulsion systems [3], leading to loss of properties and even lives.

Figure  1.  Examples of pipeline damage caused by hydraulic shock (Source: https://en.wikipedia.org/wiki/Water hammer)

DownLoad: Full-Size Img PowerPoint

Numerous methods for solving hydraulic shock model which consists of hyperbolic PDEs can be divided into three categories: graphical methods, analytical methods, and numerical methods. The former two methods require many simplifying assumptions and mostly focus on the linear pipeline transmission PDE model. They especially cannot deal with the cavitation caused by negative pressure [4]. Although these two methods are very straightforward and simple, their applications are greatly limited in the real world. Numerical methods for simulating the hydraulic shock process mainly include methods of characteristics, finite volume methods, wave characteristics methods, impulse response methods, fluid-structure interaction methods [5], and even some computational packages (e.g., COMSOL Multiphysics, FINICS, etc.). In this paper, the method of lines is used for the reduction the hydraulic shock PDEs to time varying ODEs, which is easier to understand and implement compared with other methods. Then, the ODE system is approximated through the orthogonal collocation method to obtain a spatial-temporal discretization model (larger algebraic equations). Finally, the original PDE-constrained optimization problem reduces to a nonlinear programming problem, which can be solved with many existing high-quality optimization solvers, such as sequential quadratic programming (SQP).

To protect a pipeline system from hydraulic shock effects, various protection strategies can be employed, including passive protection and active control strategies. For passive protection strategies, specialized devices are set up to enhance the ability to resist pressure shock or reduce changes in the flow rate as much as possible. For instance, material and thickness of the pipeline that can be altered to withstand high pressure shocks [6]. Pipeline networks are designed by appropriate procedures (e.g., arranging the larger pipelines in loops that supply shorter and smaller run-out pipeline branches). Water towers (generally used in many drinking water systems) help maintain steady flow rates and trap large pressure fluctuations. Other pipeline equipments, such as surge tanks, air chamber, and non-return valves, are installed for precautions against pressure shock [7]. However, these strategies described above strongly depend on the specific details of each pipeline system and on the experience of the designer. Once the design of the pipeline system and structure is complete, it is still incapable of facing a variety of conditions during the operation of the pipeline due to lack of flexibility. Hence, advanced control strategies may be used as an alternative or as a supplement to various passive protection measures.

Mitigation of hydraulic shock can be considered as a PDE control problem. Generally, there are two categories of active control strategies for dealing with this problem: boundary stabilization control and optimal control [8], [9]. The boundary stabilization control can be divided into two streams of approaches based on their characteristic method, including the Lyapunov function method [10] and the backstepping technique [11]. For the Lyapunov function method, a boundary feedback control design for open canal networks is proposed in [12], where the stabilizing boundary control law is developed by linearizing the boundary conditions. A Lyapunov function is constructed based on the weighted sum of entropies for the open channels in cascade [13]. Global stability of a version of the rotating Couette flow is examined by the sum-of-squares of polynomials approach [14]. The boundary feedback stabilization for hydraulic shock suppression is studied in [15], which can guarantee the asymptotic stability towards the steady state. For the backstepping technique, boundary stabilization for a quasilinear $2 \times 2$ system of first-order hyperbolic PDEs is considered in [16], which achieves H$^2$ exponential stability of the closed loop system. Exponential stabilization of the linearized Saint-Venant-Exner model of water dynamics in a sediment-filled canal with a single boundary actuation located at the downstream gate is achieved in [17]. Moreover, the stability of receding horizon boundary controls for scalar nonlinear hyperbolic systems is established in [18]. The control of oscillating modes modelled by the linearized Saint-Venant equations in open-channels is investigated in [19]. The works of boundary stabilization control described above mostly focus on quasilinear hyperbolic PDEs. In this paper, we consider a computational optimal control of the nonlinear hydraulic shock model instead of a feedback stabilizing controller.

Since the basic characteristics of pressure and fluid flow in the pipeline are nonlinear, nonuniform, and coupled with each other, it is difficult to establish a relationship between the valve action and the responses of the hydrodynamic variables from an optimal control point of view. An analytical expression for the rate of change of relative water head in the S-shaped region of a pumped-storage station is deduced based on a simplified mathematical model followed by a two-phase guide-vane closing scheme are proposed in [20]. Several valve closure strategies based on the known valve performance curve with unknown valve closure curve are explored in [21]. Hybrid logistic optimization of pump locations and operation based on the rainfall-runoff model and the nonlinear pump model in urban drainage systems are described in [22]. The slow unsteady-incompressible hydraulics of pressurized pipe networks are simulated and the corresponding controlled operations of the turning pumps on and off are discussed in [23]. With the advent of modern control theories and numerical methodologies, large nonlinear optimization has made the control solution of nonlinear pipeline fluid flow possible. An effective approach based on the dynamic optimization strategy can be applied to implement optimal boundary control of flows in pipelines. The control vector parameterization (CVP) method [24]–[27] for hydraulic shock suppression based on the optimize-then-discretize (OTD) method [28] and discretize-then-optimize (DTO) method [24], [29] is proposed. In the framework of the OTD method for CVP, gradient formula and the costate PDEs with optimality conditions are derived directly via the variational method based on original PDEs. However, the convergence of the OTD method for CVP is blank. Moreover, the gradients of the time-scaling terms are very complicated and computational accuracy is hard to guarantee. In the framework of the DTO method for CVP, although the convergence analysis for ODEs is complete, the CVP method and time-scaling method leads to the iterative computation cost for solving the larger scale sensitivity equations (ODE equations), which is time consuming. In practice, we should design the control strategy as soon as possible when facing hydraulic shock. Thus, in order to speed up optimization efficiency and maintain certain calculation accuracy, we focus on optimal boundary control for pressure suppression in the pipeline system via a full parameterization (FP) method [30], which may have advantages in real-time. To our best knowledge, this work is seldom studied, and offers great promise in pressure suppression of fluid transport pipelines.

The main contributions of this paper are summarized as follows: 1) A new PDE-constrained optimization for pressure suppression in fluid transport pipelines is studied. The results can further explore the applications in the pipeline network scheduling. 2) A FP control strategy based on the methods of lines and orthogonal collocation is applied to pressure suppression control. 3) Numerical simulations are carried out to validate that the FP method eliminates pressure shock effectively and costs less computation time compared with the control vector parameterization (CVP) method.

This paper is organized as follows. Section II introduces a hyperbolic PDE system for fluid flow dynamics and a valve closing model for actual valve action, and then proposes an optimal control problem for pressure suppression in the pipeline during the valve closure. Through the method of lines and the orthogonal collocation for spatial-temporal discretization, the PDE optimization problem for pressure suppression is reduced to an NLP problem in Section III. Then, a gradient-based optimization framework for solving the approximate problem, which can be solved through intrinsic subroutine FMINCON on MATLAB, is proposed. Section IV shows the effectiveness of the proposed method for pressure suppression from simulation results. Furthermore, we give the comparison between CVP method and FP method for solving this optimal control problem and find that the FP method costs less computation time. Finally, Section V concludes the paper with summary comments and suggestions for future research.

II. Problem formulation

A Mathematical Model of Pipeline Fluid Flow

Consider the situation shown in Fig. 2, where a pipeline of length $L$ is used to transport fluid from a reservoir to a terminus. Through a simplified process, where we neglect the effects of viscosity, turbulence, and temperature variation, the dynamic equations for the flow rate and pressure along the pipeline can be expressed by the following hyperbolic PDE system [31], which consists of a momentum equation and a continuity equation

Figure  2.  General layout of the pipeline system described in Section II-A

DownLoad: Full-Size Img PowerPoint

subject to the boundary conditions

where $l\in[0,L]$ stands for the spatial variable; $t\in[0,T]$ stands for the time variable; $q(l,t)$ denotes the flow rate; $p(l,t)$ denotes the edge pressure; $D$ denotes the diameter of the pipeline; $c$ denotes the wave velocity; $f$ denotes the Darcy-Weisbach friction factor; $S$ denotes the cross sectional area and $\rho$ denotes the flow density. Moreover, $P$ is the pressure generated by the reservoir (a given constant) and $u(t)$ is the terminal flow rate that represents the boundary control.

Initially, the pipeline is in steady state during the operation process. Thus, the initial flow rate of the system (1) is

where $Q$ is the steady flow rate (a constant). Substituting this condition into (1a), we can obtain

Moreover, a boundary control is required to satisfy the following bound constraint

where $U_{\max}$ denotes the maximum flow rate. Since we require the valve to be completely closed at the terminal time,

Our interest lies in modeling the fluid flow along the pipeline during the valve closure period.

B Valve Closing Model

It is not the ultimate goal to get the curve of the terminal flow rate $u(t)$ during the valve closure period since the direct control implementation is the valve closure in the actual operation. Thus, the inversion of the valve closure process, or the relationships between the valve relative opening $\alpha(t)$ and time $t$, should be demonstrated.

We denote $S(t)$ as the cross section during the valve closure period. Thus, the terminal flow rate $u(t)$ (the boundary control) is described as [32]

where $C(t)$ is the valve flow coefficient. When the valve is fully open, the valve flow coefficient $C(t)$ and the cross section $S(t)$ are expressed as $C_d$ and $S$. Substituting those into (7), we obtain

Based on (3) and (4), we divide (7) by (8) and obtain

Generally, the valve manufacturers provide the relationships between ${C(t)}/{C_d}$ and $\alpha(t)$, where $\alpha(t) = {S(t)}/{S}$. For example, in [33], the relationships between the butterfly valve relative opening $\alpha(t)$ and the valve relative flow coefficient ${C(t)}/{C_d}$ are shown in Table I.

Table  I.  The Relationships Between the Butterfly Valve Relative Opening and the Valve Relative Flow Coefficient

$\theta$	$\alpha(t)$	${C(t)}/{C_d}$
0	1	1
10	0.820	0.899
20	0.645	0.682
30	0.500	0.459
40	0.390	0.288
50	0.295	0.175
60	0.200	0.100
70	0.125	0.056
80	0.060	0.028
90	0.001	0.006

 | Show Table

DownLoad: CSV

Since the relationships between $\alpha(t)$ and ${C(t)}/{C_d}$ can be obtained, the valve closure operation depends on the boundary control $u(t)$ and corresponding pressure $p(L,t)$. If we obtain the boundary control $u(t)$, the valve relative opening $\alpha(t)$ can be obtained.

C The Optimal Boundary Control Problem

Since closing the valve suddenly can cause severe hydraulic shock effects, the boundary control $u(t)$ must be manipulated carefully to minimize such effects. Thus, we consider the following objective function, which was also used in references [34], [35]:

where $\bar{P}$ is a given constant and $P$ is the target pressure profile along the pipeline, since the final state of the pressure along the pipeline, is the same as the pressure generated by the reservoir. Thus, our optimal boundary control problem is defined as follows.

Problem P₀: Given the system (1a), (1b) with boundary conditions (2) and initial conditions (3), (4), choose the boundary control $u(t):[0,T] \to \mathbb{R}$ to minimize the objective function (10) subject to the bound constraint (5) and the terminal control constraint (6).

III. Spatial-temporal Discretization

A The Method of Lines

In order to transfer the pipeline transmission PDE model into algebraic equations, we first use the method of lines to approximate the PDE model by an ODE model [36]. The total pipeline is divided into $N_S$ equally-spaced subintervals $\left[ {l_{i-1}, l_{i} } \right], i = 1,\dots ,N_S$, where $l_0 = 0$ and $l_{N_S} = L$. Let

The $l_0, l_1, ..., l_{N_s}$ along the pipeline are shown in Fig. 3.

Figure  3.  Pipeline spatial discretization using the method of lines

DownLoad: Full-Size Img PowerPoint

Based on the definitions of $q^i(t)$ and $p^i(t)$, we can obtain the following finite difference approximations:

where $\Delta l = L/N_S$. Substituting these approximations (11a) and (11b) into the pipeline transmission PDE model (1a) and (1b) yields

By virtue of the definitions of $q^i(t)$ and $p^i(t)$, boundary conditions (2) become

And (3), (4) can be written as

We apply the Simpson’s rule [37] to approximate (10) and obtain

Since the objective function contains an integral form, we transfer it into a Mayer form by introducing a new state variable $z(t)$, where $z(0) = 0$ and $z(T) = J$

B Orthogonal Collocation Method

We approximate the ODEs through orthogonal collocation method, where those ODEs can be solved at selected points in the time domain. This method discretizes not only boundary control $u(t)$ but $q^i(t), i = 0,...,N_S-1$, $p^i(t), i = 1,...,N_S$, $z(t)$ as well, to obtain algebraic equations. This process is equivalent to a fully implicit Runge-Kutta method with excellent stability and high-order accuracy [30].

First, the total time domain $[0,T]$ is divided into $N_T$ time subintervals to obtain

where $m = 1,\dots ,N_T$ and $h_m$ is the length of $m$th time subinterval, which should be optimized. In each time subinterval, we take $(K_S+1)$ collocation points as state decision variables. Moreover, since the boundary control $u(t)$ (terminal flow rate) should be continuous, we take $(K_U+1)$ collocation points as control decision variables, where $K_U\leq K_S$. Thus, the state variables $q^i(t), i = 0,...,N_S-1$, $p^i(t), i = 1,...,N_S$, $z(t)$ can be expressed as follows utilizing piecewise Lagrange polynomials

where $t\in [t_{m},t_{m+1}]$, $\tau \in[0,1]$, $m = 1,2,\ldots ,N_T,$ $q^i_{mn} = q^i(t_{mn}),$ $p^i_{mn} = p^i(t_{mn}), z_{mn} = z(t_{mn})$ and $t_{mn} = t_{m}+h_{m}\tau_n$. The corresponding control $u(t)$ becomes

where $t\in [t_{m},t_{m+1}]$, $\tau \in[0,1]$, $m = 1,2,...,N_T,$ $u_{mj} = u(t_{mj})$ and $t_{mj} = t_{m}+h_{m}\tau_j$. The collocation points for state variable $z(t)$ are shown in Fig. 4 ($K_S = 2$). The collocation points for $u(t)$ is also shown in Fig. 4 ($K_U = 2$). Generally, Gauss-Legendre points $\tau_k, k = 1,...,K_S$ for state variables are chosen to satisfy orthogonal properties, where $\tau_0 = 0$ and $\tau_{K_S+1} = 1$. Also, $\tau_j, j = 1,...,K_U$ for control variables are chosen to satisfy orthogonal properties, where $\tau_0 = 0$ and $\tau_{K_U+1} = 1$. Moreover, since state variables (pressure, flow rate and objective function) and boundary control variables should be continuous, we can obtain the following equality constraints

Figure  4.  Collocation points for the state variable $z(t)$ and the control variable $u(t)$ ($K_S = 2, K_U = 2$)

DownLoad: Full-Size Img PowerPoint

where $m = 1,2,...,N_T-1$. The time derivative of the state variables can be represented as a Lagrange polynomial with $K_S$ collocation points. Based on the Lagrange polynomial (18) (19) and boundary condition (13), the ODEs (12) and (16) become

where

and

where

and

The initial conditions (14) become

The control constraint (5) should be bounded at each time subintervals. Thus, it can be written as

and the terminal control constraint (6) changes to

Moreover, the time subintervals should be satisfied with the following equality constraint

Finally, the objective function becomes

Remark 1: When $K_U = K_S$, the KKT conditions of the FP method are consistent with the optimality conditions of the discretized variational problem, and convergence can be guaranteed. For more details on the proof, please refer to [30]. For $K_U<K_S$, the KKT conditions of the FP method are not consistent with those of the discretized variational problem, thus we will study its convergence in the future.

Problem P_N: Given the algebraic system (21) and (23) with initial conditions (25), choose a control parameter vector $u_{mj}, m = 1,2,...,N_T, j = 1,...,K_U,$ with bound constraints (26) and terminal control constraints (27), and time subinterval variables $h_m, m = 1,2,...,N_T$ with equality constraint (28) such that the objective function (29) is minimized.

C Solving the Approximate Problem $P _N$

The approximate problem P$_N$ defined above is a nonlinear optimization problem in which a finite number of decision variables $u_{mj}, m = 1,2,...,N_T, j = 1,...,K_U$ and $h_m, m = 1,2,...,$$K_T$, need to be chosen to minimize the objective function subject to a set of constraints. For this approximate problem, the objective function and constraints are explicit functions of the decision vector. Thus, Problem P$_N$ can be now solved by serval gradient-based optimizers, such as interior-point method, SQP, active-set method, and some mature optimization softwares (e.g., IPOPT, SNOPT, MINOS, etc.) which can handle the NLP problem with intricate equality constraints, inequality constraints, and gradients of those.

We give the main steps for solving the approximate problem P$_N$. First, the total pipeline and the total time are divided into $N_S$ equally-spaced subintervals and $N_T$ time subintervals, respectively. The parameters $N_S$, $N_T$, $K_S$ and $K_U$ are set and the state variables and control variable are discretized as (18) and (19). The initial conditions are given as (25). Second, we choose the initial guess for the control parameters $u_{mj}, m = 1,2,...,N_T, j = 1,...,K_U,$ and time subinterval variables $h_m, m = 1,2,...,N_T$. Then, larger algebraic equations are constituted in (21) and (23). Finally, the objective and its gradient corresponding to the control parameters and time subinterval variables are obtained. The gradient-based optimizer performs an optimality test: if $u_{mj}$ and $h_m$ are optimal, then stop; otherwise, calculate a search direction and perform a line search to update the control vector and go to Step 2. Note that Step 3 can be implemented automatically by the gradient-based optimizer. Moreover, A gradient-based optimization framework for solving Problem P$_N$ is shown in Fig. 5.

Figure  5.  A gradient-based optimization framework for solving Problem P$_N$

DownLoad: Full-Size Img PowerPoint

IV. Numerical Simulations

For the simulation, FP method with piecewise-linear control is applied for pressure suppression in fluid transport pipelines. The pipeline parameters are given as: the total pipeline length $L = 1000 \;{\rm m}$, the diameter $D = 100\; {\rm mm}$, the Darcy-Weisbach friction factor $f = 0.03$, wave velocity $c = 1200 \;{\rm m/s}$, and the density of the liquid $\rho = 1000\;{\rm kg/m^3}$. As for the initial state of the pipeline, we can define the fluid rate $q = 1.57 \times 10^{-2} \;{\rm m/s}$ and the head of the pressure $P = 2 \times 10^6$ Pa.

We set $N_S = \{8,10,12\}$ for the spatial discretization of the pipeline and $N_T = 10$ for the temporal discretization. Using the FP method, we choose $K_U = 1$ for the piecewise-linear control and $K_S = 3$. Regarding the collocation points for state variables and control variable, we choose the Gauss–Legendre points $\tau_1 = 0.1127$, $\tau_2 = 0.5000$, and $\tau_3 = 0.8873$ for state variables, and $\tau_1 = 0.5000$ for the control variable. Moreover, the terminal time $T$ is set as 10 s, a constant $\bar{P}$ is set as $1 \times 10^6\;{\rm Pa}$.

Our numerical simulation study was carried out within the MATLAB programming environment running on a personal computer with the following configuration: Intel (R) Core (TM) i5-2320 CPU, 3.00 GHz CPU, 4.00 GB RAM, 64-bit Windows 7 Operating System, MATLAB version is R2010b. Our MATLAB code implements the gradient-based optimization procedure in Fig. 5 via the intrinsic subroutine FMINCON of MATLAB.

The initial guess for the boundary control is chosen as $-1.57 \times 10^{-2}$, which is consistent with the initial state of the flow rate, while the initial guess for the time variables is chosen as $h_m = 1, m = 1,2,\ldots,10$ which denote equidistant time subintervals. The objective function value of the initial guess, which refers to the constant closure-rate strategy, is $8.5461\times 10^{-3}$. Furthermore, we apply the CVP proposed in [29] for comparison and set the same initial guess with the FP method. Table II gives the CPU time and objective values of the FP method and the CVP method under $N_S = \{8,10,12\}$ for the spatial discretization. It shows that the computation time of the FP method is less than that of the CVP method, while the objective function of the latter is slightly less. Thus, use of the FP method is advantageous in real-time. Since the CVP method should compute $(2N_S+1) \times (2N_T+1)$ ODE equations (state and sensitive equations) with $2N_T$ decision variables in each optimization iteration to obtain the precise gradient information, computation time is very high. Meanwhile, FP method only requires solving $(2N_S+1) \times N_T\times K_S$ algebraic equations with $2N_T$ decision variables instead. For the less objective values of the CVP method, it may be more precise to use gradient information in the temporal domain, which is solved by the Runge-Kutta method instead of the Lagrange polynomial. Moreover, the objective functions of the FP and CVP methods are far less than those of constant closure-rate methods. The optimal boundary control curve ($N_S = 12$) of the FP method is presented in Fig. 6, which satisfies the control terminal constraint. Although control curves of the FP and the constant closure-rate methods are close, the corresponding objective function values are different. The computed optimal time subinterval values of the FP method are shown in Table III. The optimized time subinterval values adaptively select the best switching time instants, which result in a better control approximation. Fig. 7 shows the pressure evolution along the pipeline corresponding to the FP method ($N_S = 12$). It can be seen that, when the boundary is implemented in Fig. 6, the spatial-temporal pressure distribution fluctuates around the target pressure profile. Due to fluid physics in the pipeline, the pressure shock and oscillation can not be eliminated completely at terminal time $T$. Finally, when we have obtained the optimal boundary control $u(t)$ and corresponding pressure profile $p(L,t)$, the valve relative opening $\alpha(t)$ can be obtained by (9).

Table  II.  CPU Time and Objective Value of the FP Method and the CVP Method Under $N_S = \{8,10,12\}$ for the Spatial Discretization

Method	NS	CPU time (s)	Objective value
FP	8	396	4.829 $\times 10^{-3 }$
CVP		624	4.694 $\times 10^{-3 }$
FP	10	547	4.772 $\times 10^{-3 }$
CVP		954	4.571 $\times 10^{-3 }$
FP	12	768	4.643 $\times 10^{-3}$
CVP		1459	4.458 $\times 10^{-3}$

 | Show Table

DownLoad: CSV

Table  III.  Optimal Time Subinterval Parameters of the FP Method ($N_S = 12$)

$m$	$h_m$	$m$	$h_m$
1	0.7031	6	0.6313
2	0.6995	7	1.9160
3	1.2866	8	1.0457
4	0.6484	9	1.1537
5	1.3853	10	0.5301

 | Show Table

DownLoad: CSV

Figure  6.  Optimal boundary control curve corresponding to the FP method and constant closure-rate method ($N_S = 12$)

DownLoad: Full-Size Img PowerPoint

Figure  7.  Pressure evolution corresponding to the FP method ($N_S = 12$)

DownLoad: Full-Size Img PowerPoint

In the case study, we use the intrinsic subroutine FMINCON of MATLAB to complete FP method on the optimal boundary control problem for pressure suppression since these two methods can be performed in this software platform. Readers may wonder whether the computation time is long. There are several methods to improve the computation efficiency. For example, we have transferred the optimization problem of the pipeline PDE system for pressure suppression into an NLP problem. This NLP problem could be solved by some specialized softwares very effectively, such as AMPL, GAMS combined with optimizer tool (e.g., IPOPT, SNOPT), which can greatly reduce the computation time.

Although this NLP problem can be effectively solved by certain specialized software, the computation is not cheap enough for real-world applications, which may not meet requirements of real-time optimization. Thus, a real-time implementation that takes less computation and execution time should be designed. A more effective method to improve the real-time performance is to built databases which record many optimal solutions under different physical parameters. Solution to the closest recorded case is extracted from the database and utilized as the on-site initial guess in practical applications. We test the computation cost under a suboptimal initial guess, which shows the optimal result can be obtained in around 9s on MATLAB. Of course, this result is not comprehensive and sufficient. Moreover, inspired by the operation of the beam pumping units [38], [39], the aim of this paper is to apply the proposed method into practice, we will study in the future.

V. Conclusion and Future Work

This paper has presented an FP method for solving the PDE optimal control problem during valve closure periods towards pressure suppression. This method is based on a combination of the spatial discretization use the method of lines and time discretization by the orthogonal collocation method. The simulation results demonstrate the FP method is superior to the CVP method in terms of computation time. The simulation results demonstrate that the proposed method is effective toward pressure suppression. Since real-time implementation is more practical in the real-work by suitable scheduling of the valve operation sequences, the databases which record many optimal solutions under different physical parameters should be effectively built. Moreover, an improved predictive control strategy, dynamic programming with guaranteed satisfaction of path constraints, and a multi-objective optimization method for the nonlinear dynamic optimization are considered in [40]–[42], which can be extended to the spatial-temporal hydraulic shock model. These directions will be explored in future work.