# Discussion on Muskingum versus Integrator-Delay Models for Control Objectives

Authors: Yolanda Bolea; Vicenç Puig; Antoni Grau
Journal: Journal of Applied Mathematics
 (2014)

---

## Abstract

A comparative study about two models, Muskingum and integrator-delay (ID) models, for canal control is presented. The former is a simplified hydrological model which is very simple and extensively used in hydraulic engineering for simulation and prediction. The latter is also a model with physical meaning and is widely used for irrigation canals control. Due to a lack of general awareness of Muskingum prediction model in regulation from the control community, authors present this comparative study with the ID control model. Both models have been studied and analyzed for control purposes. This study has been carried out and validated in a real irrigation canal, at Aghili irrigation district in Iran, using two traditional control approaches, PID with feedback and predictive control. The results demonstrate the advantages and drawbacks of both models, showing the benefits and limitations of using the widespread Muskingum model among the hydraulics scientific community for control design.


---

## Body

## 1. Introduction

Management of open-flow canal systems requires accurate control models of flow transfer. Open-flow canals are large parameter-distributed systems that can be described with Saint-Venant equations [1, 2]. These nonlinear partial differential equations (PDE) represent water dynamics in a precise and complete manner and, for an arbitrary geometry, there is no analytical solution. Usually numerical methods [3, 4] have been used to obtain a solution of the Saint-Venant equations. The more typical numerical methods are the characteristic method [5], the finite difference method [6], and Preissman implicit scheme [7–9]. The implicit numerical schemes are simpler than the other approaches. But the stability of their solution cannot be guaranteed and depends on the discretization time. Taking into account this fact, usually Muskingum model has been used by hydraulic engineers as a prediction model in rivers and in irrigation canals [10, 11] but so far not as a control model. Nevertheless, Saint-Venant model is not useful for automatic control purposes due to its complexity (high number of states and high computational load, etc.). Unfortunately, in control design or in the computation of control actions, the complexity of the characteristic model is directly proportional to the complexity of the control techniques and their implementation. Modeling these systems for control design is therefore not so easy to be devised, although it is a crucial step. A classical way to control irrigation canals is to design controllers based on linear models.Therefore, in the literature several linear control models have been already proposed for open-flow canal systems. Corriga et al. [12, 13] proposed a model for open canal networks related to the elevation of the gate. Papageorgiou and Messmer [14, 15] and Ermolin [16] proposed an approximation for uniform flow canals derived from linearized Saint-Venant equations. Malaterre [7] used a state-space model for canal systems derived from discretized Saint-Venant equations. Baume et al. [17] obtained an infinite dimensional state-space model from linearized Saint-Venant equations. When the canal has prismatic geometry and it is in permanent regime, Baume’s model has the advantage that the solution of the linearized Saint-Venant equations is exact. But usually, real canals do not fulfill the required conditions assumed in the mentioned works. To palliate these problems, Schuurmans et al. [18–20] formulated an analytical model denoted as integrator delay (ID) that is composed of two parameters: an integrator and a delay. This model is an approximation that relates input and output flow including backwater effects when a reservoir-delay model is used. It is very popular among the control community for its easiness describing the more essential aspects of the canal frequential behavior. However, in [21], it is demonstrated that the ID model does not represent accurately the dynamics of a canal for medium and high frequencies. Another extended model for control purposes is the Hayami model [22, 23]. It is based on a first order linearization of the diffusive wave equation around the reference flow. This model is represented by a second order transfer function with delay structure. The model behaves accurately for abrupt and unexpected changes of the reference signal. Its main drawback is its usefulness only for canals with a specific geometry such as those with a reduced set of diffusion and celerity coefficients and canal length.To overcome the drawbacks of the mentioned models so far, some relatively recent works have been developed. Improvements to the ID model have been also proposed in [24] and its LPV extension [25]. Moreover, in [26, 27], an IDZ (integrator delay zero) model is proposed. This novel model extends the ID model by adding a zero in the high frequencies leading to a better fit in this range. This fact improves the accuracy in time domain simulations, and all model parameters can be analytically computed leading to a simpler implementation. This modelling approach has been also used to generate state-space MIMO models in [28–31]. In [32], a new computational method is provided in order to obtain a frequency domain model of Saint-Venant equations linearized around any stationary regime. In [33], a mathematical model of a gate to automatically control the upstream water level is proposed in order to operate under free and submerged flows conditions (called the Vlugter gate). Munier and his colleagues, in [34], proposed a new model, called LBLR for Linear Backwater Lag-and-Route, which approximates the Saint-Venant equations linearized around a nonuniform flow in a finite channel (with a downstream boundary condition) taking into account the backwater effect. Next in [35], a new approach to compute the response time is proposed, accounting explicitly for the backwater and the feedback effects due to the downstream cross structure. The method provides a distributed analytical expression of the response time as a function of the characteristics of the canal (geometry and roughness) and the downstream cross structure.As an alternative to the above approaches, in [36–39], linear black-box models for the canal control design have been used. These models are estimated by classical linear estimation methods [40]. They do not have hydraulic meaning and they are based only on experiments. In this case, the main drawback is the necessity of a high number of experiments, especially in the multivariable case. For this reason, some authors carried out a subspace identification [41, 42].Interestingly, there is a model widely used in hydraulic engineering for prediction that has been barely used as a control model for controller design. This is the Muskingum model. After his experiments in the Muskingum River, McCarthy proposed a hydrological model suitable only in rivers (Muskingum model, [10, 43]) operating with input and output flows. The model does not require a geometrical definition for the riverbed because the geometry is included in its parameters.Muskingum model is a hydrological model and is widely used for simulation by hydraulic engineers due to its simplicity and its good results for prediction; nevertheless, it has been almost never applied in the control of flow-open channels. On the other hand, ID model is a mere control model widely used in this kind of systems and with satisfactory results. Both models are simple and simplified and, moreover, they have physical meaning. And this is the motivation that pushed this research to analyze and validate this hydrological model in the canal control and compare it with ID model. Muskingum model is very interesting for control because it is a discrete model that can be used with the most known real canal control methods based on PID and predictive control. In engineering, Muskingum model is not used mainly because it is an unknown approach and the objective of this work is also to reveal, investigate, and study the advantages and drawbacks of Muskingum model (typically used in prediction) as an open-flow canal control model. To reach this objective, a comparative study between both models (ID model and Muskingum model) has been carried out. It is well worth to note that in [44], the Muskingum model is used in predictive control. It would be of interest that hydraulic engineers as well as control engineers could use the same mathematical tools to share experiences and knowledge, because at the end they are performing complementary works in water resource management. Moreover, both communities could approach and create synergies with a higher efficacy and efficiency.The main contribution of this research is to analyze and study whether Muskingum hydrological model is suitable as control model for typical controller design in water management as feedback control with PID and predictive control. The present paper is organized as follows. Section2 gives a system description of Aghili irrigation canal used in our study. In Section 3, a brief explanation of Muskingum model and its relationship with the ID model is presented. Specifically, the Muskingum prediction model has been studied as a control model. To do so, it has been analyzed in open loop (in continuous time as well as in discrete time) and in closed loop (considering the effect of ZOH and sample time). In Section 4, a comparison is also performed with the ID model widely used in controller design. In Section 5, the closed-loop behaviour in continuous and discrete time is also analyzed specially when using such a model to design a feedback controller. In Section 6, both models are tested and validated in Aghili irrigation canal using two control methods: PID control with a SP scheme and predictive control. Finally, the conclusions of using the Muskingum and ID models for control are presented.


## 2. Description of the Irrigation Canal

In this section, a real irrigation canal system (Aghili canal) is described. The complete behaviour of the water in this system is represented by Saint-Venant equations as a simulation model.


### 2.1. System Description

Aghili irrigation district (AID) is located in southwest Iran, in the north of Khuzestan Province (see Figure1(a)). AID is a part of the Gotvand irrigation network. The annual (maximum versus minimum) air temperatures range from 53° to 3°C and precipitation rates from 582 mm to 152 mm, respectively. The net cultivated area in AID is about 4000 ha. The annual mean distributed water in the irrigation area is about 150 MCM (million cubic meters). AID includes a short main canal, 2 km long, along with two subsidiary canals, right branch at a length of 14.9 km and left branch at a length of 18.6 km. In this study, the 2 km canal is treated. In Figure 1(b), an aerial picture of the spillway can be seen where both branches depart from the canal. The right branch (west branch) has been studied in [45]. The total irrigation canal is composed of a main small pool equipped with an upstream sluice gate that takes water from Karun River (Figure 1(c)) with two subsidiary reaches (Figure 2(a)), as previously cited. Upstream of this gate (u) is a dam (Figure 1(d)) of constant level H
=
3.5 m and the total length of the pool (L) is 2 km (Figure 2(b)). For more information about Karun River, see [46]. This main pool is used for the control study of this paper.(a) Map of Khuzestan Province in Iran where the Aghili irrigation map is located; (b) Aghili irrigation canal and its spillway, West and East branches start at this point (32.169147 and 48.734554 coordinates) (from Google Maps); (c) aerial picture of the start of the canal at Karun River (from Google Maps); (d) sluice gate at Karun river where the 2 km canal starts.
(a)
(b)
(c)
(d)(a) Canal scheme and (b) scheme of the main reach.
(a)
(b)To reproduce the real behaviour of this irrigation canal, a simulator developed by the group of “Modelling and Control of Hydraulic Systems” at the UPC is used [47]. It solves numerically Saint-Venant equations [1, 3, 48] which accurately describe water dynamics.


### 2.2. Saint-Venant Model

Saint-Venant model can be expressed by the conservation of mass and momentum principles equations in a one-dimensional free surface flow.Conservation mass equation:(1)
∂
q
∂
x
+
∂
A
∂
t
=
0
.Momentum equation:(2)
∂
q
∂
t
+
∂
∂
x
(
q
2
A
)
+
g
A
∂
y
∂
x
-
g
A
(
I
0
-
I
f
)
=
0
,
where Q
=
Q
(
x
,
t
) is the flow (m3/s), y
(
x
,
t
) is the level (m), A
=
A
(
x
,
t
) is the cross-sectional area (m2), t is the time variable (s), x is the spatial variable (m) measured in the direction and the sense of the movement, y is the water level, g is the gravity (m/s2), I
0 is the bottom slope, and I
f is the friction slope.This pair of partial differential equations (1)-(2) constitutes a nonlinear and hyperbolic system for an arbitrary geometry and it lacks analytical solution [1, 3]. These equations can be solved by different numerical methods, such as Preissman, finite differences, and characteristics [48, 49]. The Saint-Venant model describes in detail the canal dynamics, and it is very useful in simulation. However, this model is too complex for controller design and furthermore it cannot be used for control design using analytical methods. As previously said, these equations are used as simulation model (solved numerically using the method of characteristics) to reproduce the complex dynamics of Aghili irrigation canal.


## 2.1. System Description

Aghili irrigation district (AID) is located in southwest Iran, in the north of Khuzestan Province (see Figure1(a)). AID is a part of the Gotvand irrigation network. The annual (maximum versus minimum) air temperatures range from 53° to 3°C and precipitation rates from 582 mm to 152 mm, respectively. The net cultivated area in AID is about 4000 ha. The annual mean distributed water in the irrigation area is about 150 MCM (million cubic meters). AID includes a short main canal, 2 km long, along with two subsidiary canals, right branch at a length of 14.9 km and left branch at a length of 18.6 km. In this study, the 2 km canal is treated. In Figure 1(b), an aerial picture of the spillway can be seen where both branches depart from the canal. The right branch (west branch) has been studied in [45]. The total irrigation canal is composed of a main small pool equipped with an upstream sluice gate that takes water from Karun River (Figure 1(c)) with two subsidiary reaches (Figure 2(a)), as previously cited. Upstream of this gate (u) is a dam (Figure 1(d)) of constant level H
=
3.5 m and the total length of the pool (L) is 2 km (Figure 2(b)). For more information about Karun River, see [46]. This main pool is used for the control study of this paper.(a) Map of Khuzestan Province in Iran where the Aghili irrigation map is located; (b) Aghili irrigation canal and its spillway, West and East branches start at this point (32.169147 and 48.734554 coordinates) (from Google Maps); (c) aerial picture of the start of the canal at Karun River (from Google Maps); (d) sluice gate at Karun river where the 2 km canal starts.
(a)
(b)
(c)
(d)(a) Canal scheme and (b) scheme of the main reach.
(a)
(b)To reproduce the real behaviour of this irrigation canal, a simulator developed by the group of “Modelling and Control of Hydraulic Systems” at the UPC is used [47]. It solves numerically Saint-Venant equations [1, 3, 48] which accurately describe water dynamics.


## 2.2. Saint-Venant Model

Saint-Venant model can be expressed by the conservation of mass and momentum principles equations in a one-dimensional free surface flow.Conservation mass equation:(1)
∂
q
∂
x
+
∂
A
∂
t
=
0
.Momentum equation:(2)
∂
q
∂
t
+
∂
∂
x
(
q
2
A
)
+
g
A
∂
y
∂
x
-
g
A
(
I
0
-
I
f
)
=
0
,
where Q
=
Q
(
x
,
t
) is the flow (m3/s), y
(
x
,
t
) is the level (m), A
=
A
(
x
,
t
) is the cross-sectional area (m2), t is the time variable (s), x is the spatial variable (m) measured in the direction and the sense of the movement, y is the water level, g is the gravity (m/s2), I
0 is the bottom slope, and I
f is the friction slope.This pair of partial differential equations (1)-(2) constitutes a nonlinear and hyperbolic system for an arbitrary geometry and it lacks analytical solution [1, 3]. These equations can be solved by different numerical methods, such as Preissman, finite differences, and characteristics [48, 49]. The Saint-Venant model describes in detail the canal dynamics, and it is very useful in simulation. However, this model is too complex for controller design and furthermore it cannot be used for control design using analytical methods. As previously said, these equations are used as simulation model (solved numerically using the method of characteristics) to reproduce the complex dynamics of Aghili irrigation canal.


## 3. Physical Canal Models: Muskingum Model versus ID model

### 3.1. Muskingum Model

The Muskingum model [1, 10] relates the inflow (I) with the outflow (Q) of a free surface hydraulic system. This model uses the following continuity equation and a linear storage-discharge relationship:
(3)
I
(
t
)
-
Q
(
t
)
=
d
V
(
t
)
d
t
,
(4)
V
(
t
)
=
K
(
χ
I
(
t
)
+
(
1
-
χ
)
Q
(
t
)
)
,
where V (m3) is the volume of water stored in the hydraulic system and K and χ are parameters that depend on the hydraulic system. K is considered as the average pool travel time, and χ shows the influence of the downstream condition in the hydrograph propagation with 0
≤
χ
≤
0.5 (see Figure 3).Figure 3
Prismatic storage and wedge-shaped storage in a pool.The inflowI is the input variable, also noted as q
ups; the outflow Q is the output variable, also noted q
dns. The Muskingum model works only with upstream and downstream flows and it is not possible to obtain the water level with any direct relationship. Normally, an empirical approximation is used (e.g., Manning, Chezy or Darcy-Weisbach; in this paper, the Manning approximation is used) to relate the flow and level, adding imprecision to the results. Replacing (3) with (4) and using the Laplace transform, the relation between output flow and input flow is obtained:
(5)
q
dns
(
s
)
=
G
m
(
s
)
q
ups
(
s
)
=
(
1
-
K
χ
s
)
(
K
(
1
-
χ
)
s
+
1
)
q
ups
(
s
)
.


### 3.2. ID Model

Canals are systems that can be decomposed in many dynamic elements, usually of first order, so the full model is of an order equal to the number of elements. Then, the resulting model will have an order equal to the number of pieces used to model the canal, that is, a very high order model. As discussed in [50], this very high order model would be difficult to use for control purposes but, fortunately, it is possible to approximate the behaviour of such high order processes by a system with one time constant and a dead time, that is, a FOPDT. This will be shown in the following subsection.The complete water dynamics of a single-pool irrigation canal is classically modeled with the Saint-Venant equations. However, as discussed in the previous section, for control purposes, the ID (integrator delay) model is considered as proposed by Litrico and Fromion [32] for low frequencies. The single canal reach dynamics has a relationship between downstream level y
dns and downstream flow q
dns and upstream flow q
ups for low frequencies. This relation can be approximated by
(6)
y
dns
(
s
)
=
P
1
(
s
)
q
ups
(
s
)
+
P
2
(
s
)
q
dns
(
s
)
with
(7)
P
1
(
s
)
=
e
-
τ
d
s
A
d
s
,
P
2
(
s
)
=
-
1
A
d
s
,
where τ
d is the downstream transport delay and A
d is the downstream backwater area. Taking into account the linearized relationship between the downstream flow q
dns and level y
dns in the spillway, the following equality can be established:
(8)
q
dns
(
s
)
=
λ
y
dns
(
s
)
,
where λ is a constant.Combining (6) and (8), the following first order plus time delay (FOPDT) model can be obtained:
(9)
q
dns
(
s
)
=
P
1
(
s
)
(
λ
-
P
2
(
s
)
)
q
ups
(
s
)
=
e
-
τ
d
s
(
A
d
λ
s
+
1
)
q
ups
(
s
)
=
k
(
T
s
+
1
)
e
-
τ
d
s
q
ups
(
s
)
with a gain k
=
1 and a time constant T
=
λ
A
d.The delay,τ
d, can be estimated based on physical laws by
(10)
τ
d
=
L
v
+
c
,
where L is the canal length, v is the current water velocity, and c is the current celerity: c
=
g
y. For more details about this modelling approach, see [22–24, 26, 32].


### 3.3. Relation between Both Models

Comparing the ID and Muskingum model structure, it can be shown that Muskingum model can be obtained from the ID model through an approximation using the McLaren approach [51]. Then, if G
m is the open-flow canal transfer function, we obtain
(11)
G
m
(
s
)
=
1
T
s
+
1
e
-
τ
s
≈
1
-
τ
s
T
s
+
1
.
Although the delay is neither close to zero nor it can be neglected considering the system dynamics, the error produced by the approximation is not significant when G
m is discretized by a convenient sampling time; see Section 3.2. This approximated transfer function has the same structure as that of the hydrological original continuous Muskingum model, considering τ
=
K
χ and T
=
K
(
1
-
χ
).


## 3.1. Muskingum Model

The Muskingum model [1, 10] relates the inflow (I) with the outflow (Q) of a free surface hydraulic system. This model uses the following continuity equation and a linear storage-discharge relationship:
(3)
I
(
t
)
-
Q
(
t
)
=
d
V
(
t
)
d
t
,
(4)
V
(
t
)
=
K
(
χ
I
(
t
)
+
(
1
-
χ
)
Q
(
t
)
)
,
where V (m3) is the volume of water stored in the hydraulic system and K and χ are parameters that depend on the hydraulic system. K is considered as the average pool travel time, and χ shows the influence of the downstream condition in the hydrograph propagation with 0
≤
χ
≤
0.5 (see Figure 3).Figure 3
Prismatic storage and wedge-shaped storage in a pool.The inflowI is the input variable, also noted as q
ups; the outflow Q is the output variable, also noted q
dns. The Muskingum model works only with upstream and downstream flows and it is not possible to obtain the water level with any direct relationship. Normally, an empirical approximation is used (e.g., Manning, Chezy or Darcy-Weisbach; in this paper, the Manning approximation is used) to relate the flow and level, adding imprecision to the results. Replacing (3) with (4) and using the Laplace transform, the relation between output flow and input flow is obtained:
(5)
q
dns
(
s
)
=
G
m
(
s
)
q
ups
(
s
)
=
(
1
-
K
χ
s
)
(
K
(
1
-
χ
)
s
+
1
)
q
ups
(
s
)
.


## 3.2. ID Model

Canals are systems that can be decomposed in many dynamic elements, usually of first order, so the full model is of an order equal to the number of elements. Then, the resulting model will have an order equal to the number of pieces used to model the canal, that is, a very high order model. As discussed in [50], this very high order model would be difficult to use for control purposes but, fortunately, it is possible to approximate the behaviour of such high order processes by a system with one time constant and a dead time, that is, a FOPDT. This will be shown in the following subsection.The complete water dynamics of a single-pool irrigation canal is classically modeled with the Saint-Venant equations. However, as discussed in the previous section, for control purposes, the ID (integrator delay) model is considered as proposed by Litrico and Fromion [32] for low frequencies. The single canal reach dynamics has a relationship between downstream level y
dns and downstream flow q
dns and upstream flow q
ups for low frequencies. This relation can be approximated by
(6)
y
dns
(
s
)
=
P
1
(
s
)
q
ups
(
s
)
+
P
2
(
s
)
q
dns
(
s
)
with
(7)
P
1
(
s
)
=
e
-
τ
d
s
A
d
s
,
P
2
(
s
)
=
-
1
A
d
s
,
where τ
d is the downstream transport delay and A
d is the downstream backwater area. Taking into account the linearized relationship between the downstream flow q
dns and level y
dns in the spillway, the following equality can be established:
(8)
q
dns
(
s
)
=
λ
y
dns
(
s
)
,
where λ is a constant.Combining (6) and (8), the following first order plus time delay (FOPDT) model can be obtained:
(9)
q
dns
(
s
)
=
P
1
(
s
)
(
λ
-
P
2
(
s
)
)
q
ups
(
s
)
=
e
-
τ
d
s
(
A
d
λ
s
+
1
)
q
ups
(
s
)
=
k
(
T
s
+
1
)
e
-
τ
d
s
q
ups
(
s
)
with a gain k
=
1 and a time constant T
=
λ
A
d.The delay,τ
d, can be estimated based on physical laws by
(10)
τ
d
=
L
v
+
c
,
where L is the canal length, v is the current water velocity, and c is the current celerity: c
=
g
y. For more details about this modelling approach, see [22–24, 26, 32].


## 3.3. Relation between Both Models

Comparing the ID and Muskingum model structure, it can be shown that Muskingum model can be obtained from the ID model through an approximation using the McLaren approach [51]. Then, if G
m is the open-flow canal transfer function, we obtain
(11)
G
m
(
s
)
=
1
T
s
+
1
e
-
τ
s
≈
1
-
τ
s
T
s
+
1
.
Although the delay is neither close to zero nor it can be neglected considering the system dynamics, the error produced by the approximation is not significant when G
m is discretized by a convenient sampling time; see Section 3.2. This approximated transfer function has the same structure as that of the hydrological original continuous Muskingum model, considering τ
=
K
χ and T
=
K
(
1
-
χ
).


## 4. Open-Loop Behaviour of the Muskingum Model

### 4.1. Continuous-Time Open-Loop Behaviour

The behaviour of the Muskingum model in open-loop will be compared with the ID model and with the real behaviour reproduced by means of solving the Saint-Venant equations [1] using the open-flow canal presented in Section 2 of this paper. The estimated parameters of the Muskingum model for this canal are K
=
1200 and χ
=
0.2 while in the case of ID model they are k
=
1, T
=
12.09 min, and τ
=
8.28 min.Equation (11) shows that the continuous Muskingum model is a model with a zero in the right half-plane (RHP), that is, with unstable inverse model (nonminimum phase model, NMP). This implies that for either step-input or instantaneous unit hydrograph response will present negative flows (i.e., when the input increases the output decreases) lacking physical meaning. This problem has been already reported in [52]. In Figure 4, the step-input response of the canal system for the case of the Muskingum and ID model is compared against the real behaviour simulated using the Saint-Venant equations. It can be noticed that continuous ID behaviour resembles the real one (simulated with the Saint-Venant equations). On the other hand, the nonminimum phase behaviour of the continuous Muskingum model is present at the beginning, although it reaches the correct steady state at the end.Figure 4
Real-step response (Saint-Venant equations) and estimated-step responses with continuous Muskingum model and continuous ID model.


### 4.2. Discrete-Time Open-Loop Behaviour

When obtaining the discrete-time Muskingum model to predict the behaviour of an open flow sewer in open-loop, the bilinear transform is usually used in the hydraulic literature. Then, applying the bilinear transformation with sampling timeT
s to (1) and (2), the following expression is obtained:
(12)
Q
(
k
)
=
c
0
I
(
k
)
+
c
1
I
(
k
-
1
)
+
c
2
Q
(
k
-
1
)
,
where
(13)
c
0
=
T
s
/
K
-
2
χ
2
(
1
-
χ
)
+
T
s
/
K
,
c
1
=
T
s
/
K
+
2
χ
2
(
1
-
χ
)
+
T
s
/
K
,
c
2
=
2
(
1
-
χ
)
-
T
s
/
K
2
(
1
-
χ
)
+
T
s
/
K
with the condition
(14)
2
K
χ
≤
T
s
≤
K
in order to provide hydrological meaning to the model. This selection guarantees that the sampling time T
s will be larger than the delay and the propagated hydrograph does not take negative values.Then, the following transfer function can be obtained from (12)
(15)
Q
(
z
)
I
(
z
)
=
c
0
z
+
c
1
z
-
c
2
that corresponds to a first-order linear system with a zero and a pole. According to (15), in fact the parameters will be c
0
,
c
1, and c
2 and using the expressions in (13) it will be possible to compute Muskingum parameters K and χ (bilinear transformation, see Table 1).Table 1
Muskingum parameters.
Discretization type
Muskingum coefficients
c
0
c
1
c
2
With ZOH
-
χ
(
1
-
χ
)
(
χ
1
-
χ
)
+
(
1
-
e
-
T
s
/
K
(
1
-
χ
)
)
l
l
l
l
l
l
l
e
-
T
s
/
K
(
1
-
χ
)
g
g
g
g
g
g
g
g
g
With bilinear transform
T
-
2
K
χ
2
K
(
1
-
χ
)
+
T
2
K
(
1
-
χ
)
-
T
2
K
(
1
-
χ
)
+
T
T
+
2
K
χ
2
K
(
1
-
χ
)
+
T
g
g
g
g
gIn Figure5, the open-loop response of the discrete-time Muskingum model to a unit step input is compared with the one obtained using an ID model and the real behaviour simulated using Saint-Venant equations. It can be noticed that using the sampling time suggested by rule equation (14) and the bilinear discretization method, the discrete-time Muskingum model produces good prediction results, avoiding the negative values that presented the continuous-time model.Figure 5
Real-step response (Saint-Venant equations) and estimated-step responses with continuous and discrete (bilinear,T
s
=
600 s) Muskingum models and continuous and discrete (ZOH, T
s
=
60 s) ID model.


## 4.1. Continuous-Time Open-Loop Behaviour

The behaviour of the Muskingum model in open-loop will be compared with the ID model and with the real behaviour reproduced by means of solving the Saint-Venant equations [1] using the open-flow canal presented in Section 2 of this paper. The estimated parameters of the Muskingum model for this canal are K
=
1200 and χ
=
0.2 while in the case of ID model they are k
=
1, T
=
12.09 min, and τ
=
8.28 min.Equation (11) shows that the continuous Muskingum model is a model with a zero in the right half-plane (RHP), that is, with unstable inverse model (nonminimum phase model, NMP). This implies that for either step-input or instantaneous unit hydrograph response will present negative flows (i.e., when the input increases the output decreases) lacking physical meaning. This problem has been already reported in [52]. In Figure 4, the step-input response of the canal system for the case of the Muskingum and ID model is compared against the real behaviour simulated using the Saint-Venant equations. It can be noticed that continuous ID behaviour resembles the real one (simulated with the Saint-Venant equations). On the other hand, the nonminimum phase behaviour of the continuous Muskingum model is present at the beginning, although it reaches the correct steady state at the end.Figure 4
Real-step response (Saint-Venant equations) and estimated-step responses with continuous Muskingum model and continuous ID model.


## 4.2. Discrete-Time Open-Loop Behaviour

When obtaining the discrete-time Muskingum model to predict the behaviour of an open flow sewer in open-loop, the bilinear transform is usually used in the hydraulic literature. Then, applying the bilinear transformation with sampling timeT
s to (1) and (2), the following expression is obtained:
(12)
Q
(
k
)
=
c
0
I
(
k
)
+
c
1
I
(
k
-
1
)
+
c
2
Q
(
k
-
1
)
,
where
(13)
c
0
=
T
s
/
K
-
2
χ
2
(
1
-
χ
)
+
T
s
/
K
,
c
1
=
T
s
/
K
+
2
χ
2
(
1
-
χ
)
+
T
s
/
K
,
c
2
=
2
(
1
-
χ
)
-
T
s
/
K
2
(
1
-
χ
)
+
T
s
/
K
with the condition
(14)
2
K
χ
≤
T
s
≤
K
in order to provide hydrological meaning to the model. This selection guarantees that the sampling time T
s will be larger than the delay and the propagated hydrograph does not take negative values.Then, the following transfer function can be obtained from (12)
(15)
Q
(
z
)
I
(
z
)
=
c
0
z
+
c
1
z
-
c
2
that corresponds to a first-order linear system with a zero and a pole. According to (15), in fact the parameters will be c
0
,
c
1, and c
2 and using the expressions in (13) it will be possible to compute Muskingum parameters K and χ (bilinear transformation, see Table 1).Table 1
Muskingum parameters.
Discretization type
Muskingum coefficients
c
0
c
1
c
2
With ZOH
-
χ
(
1
-
χ
)
(
χ
1
-
χ
)
+
(
1
-
e
-
T
s
/
K
(
1
-
χ
)
)
l
l
l
l
l
l
l
e
-
T
s
/
K
(
1
-
χ
)
g
g
g
g
g
g
g
g
g
With bilinear transform
T
-
2
K
χ
2
K
(
1
-
χ
)
+
T
2
K
(
1
-
χ
)
-
T
2
K
(
1
-
χ
)
+
T
T
+
2
K
χ
2
K
(
1
-
χ
)
+
T
g
g
g
g
gIn Figure5, the open-loop response of the discrete-time Muskingum model to a unit step input is compared with the one obtained using an ID model and the real behaviour simulated using Saint-Venant equations. It can be noticed that using the sampling time suggested by rule equation (14) and the bilinear discretization method, the discrete-time Muskingum model produces good prediction results, avoiding the negative values that presented the continuous-time model.Figure 5
Real-step response (Saint-Venant equations) and estimated-step responses with continuous and discrete (bilinear,T
s
=
600 s) Muskingum models and continuous and discrete (ZOH, T
s
=
60 s) ID model.


## 5. Closed-Loop Behaviour of the Muskingum Model

The use of the Muskingum model in closed-loop canal control (digital and continuous) is discussed in this section taking into account three important features for control applications such as the effect of its structure, the inclusion of a zero-order hold (ZOH) in closed loop, and the sampling time selection. This study has been carried out considering a PI controller synthesized by the pole placement method [53] and using the same canal system in case of the open-loop response (Section 4). The PI controller has been designed in order to obtain the following specifications: steady-state zero error for step responses, no overshoot, and desired time constant of 800 s.


### 5.1. Structure

As discussed in Section3.1, the Muskingum model introduces a nonminimum phase zero. The effect of such zero in closed loop is the degradation of the performance because the controller has been adjusted considering that the canal behaves as a nonminimum phase plant but the real canal does not behave in such a way. In particular, the desired closed-loop behaviour is not obtained, appearing instead as an overshoot when applied to the real canal. On the other hand, when the Muskingum model is used as the plant to simulate the behaviour of the controller that also has been designed using the Muskingum model, the closed behaviour fulfills the desired behaviour but still presents the nonminimum phase dynamics (Figure 6). As a conclusion, the effect of structure mismatch between the Muskingum model and the real canal behaviour will derive in a degraded closed-loop response including more overshoot than the specified.Figure 6
Step closed-loop response for the case I: downstream flow with PI designed by continuous Muskingum with continuous Muskingum plant; and case II: downstream flow with PI designed by continuous Muskingum with ideal ID plant.


### 5.2. Discretization Method

When using the discrete Muskingum model for designing a closed-loop digital controller, additionally to the problem of the structure mismatch, two additional issues appear; namely, the discrete equivalent Muskingum model should include the dynamics of a ZOH and second the Shannon criteria should be fulfilled when selecting the sampling time. These problems appear since, usually in hydraulics, the Muskingum model is discretized using the bilinear transformation and with a sampling time that provides hydrological meaning to the model (see Section3.1).However, when using such a model in closed-loop digital control, instead a ZOH transform should be used to take into account the dynamics of the ZOH that models the D/A conversion process. This is a standard procedure in classical digital control [54]. Then, equivalent discrete-time transfer function will have the same structure as that in the case of the bilinear transform, but the parameters of the model c
0
,
c
1, and c
2 will be different according to Table 1.If the discrete-time Muskingum model obtained by means of the bilinear transform is used to design the controller instead of the one obtained using the zero-order hold transform, the closed-loop behaviour will experiment a performance degradation when applied to a discrete-time Muskingum model using the zero-order hold transform according to the results presented in Figure7. The inclusion of the zero-order hold when applying digital control is unavoidable since the D/A converter is always present. As a conclusion, for closed-loop digital control design, the discrete-time Muskingum model should be obtained using the zero-order hold transform instead of the bilinear transform.Figure 7
Step closed-loop response for the case I: downstream flow with PI designed by continuous Muskingum with continuous Muskingum plant; case II: downstream flow with PI designed by discrete Muskingum (ZOH andT
s
=
60 s) with discrete Muskingum (ZOH and T
s
=
60 s) plant; case III: downstream flow with PI designed by discrete Muskingum (bilinear and T
s
=
600 s) with discrete Muskingum (bilinear and T
s
=
600 s) plant; case IV: downstream flow with PI designed by discrete Muskingum (bilinear and T
s
=
600 s) with discrete Muskingum (ZOH and T
s
=
600 s) plant; and case V: downstream flow with PI designed by discrete Muskingum (ZOH and T
s
=
600 s) with discrete Muskingum (ZOH and T
s
=
600 s) plant.


### 5.3. Effects of the Sampling Time

Even when discretizing the Muskingum model using the zero-order hold method, another problem appears when selecting the sampling time period satisfying the rule given in (9) in order that the discrete Muskingum model has hydrological meaning since the Shannon theorem [54] is not respected. Shannon’s theorem establishes that a good sampling time should take a value at least T
s
=
T
/
10, but according to rule (14), in the best case T
s
=
K
=
T
/
(
1
-
χ
). The selection of the sampling time following the Muskingum rule will derive again in a performance’s degradation according to Figure 7 with respect to the equivalent closed-loop continuous behaviour. In fact, in this figure it can be observed that the use of the sampling time suggested by Muskingum rule increases the overshoot. On the other hand, selecting the sampling time according to Shannon’s theorem derives in a closed-loop response that resembles that obtained when the continuous Muskingum model is used. However, even when using the ZOH transform method as a discretization procedure and the sampling time suggested by Shannon’s theorem, the problem of the structure in the Muskingum model described in Section 5.1 still remains.


## 5.1. Structure

As discussed in Section3.1, the Muskingum model introduces a nonminimum phase zero. The effect of such zero in closed loop is the degradation of the performance because the controller has been adjusted considering that the canal behaves as a nonminimum phase plant but the real canal does not behave in such a way. In particular, the desired closed-loop behaviour is not obtained, appearing instead as an overshoot when applied to the real canal. On the other hand, when the Muskingum model is used as the plant to simulate the behaviour of the controller that also has been designed using the Muskingum model, the closed behaviour fulfills the desired behaviour but still presents the nonminimum phase dynamics (Figure 6). As a conclusion, the effect of structure mismatch between the Muskingum model and the real canal behaviour will derive in a degraded closed-loop response including more overshoot than the specified.Figure 6
Step closed-loop response for the case I: downstream flow with PI designed by continuous Muskingum with continuous Muskingum plant; and case II: downstream flow with PI designed by continuous Muskingum with ideal ID plant.


## 5.2. Discretization Method

When using the discrete Muskingum model for designing a closed-loop digital controller, additionally to the problem of the structure mismatch, two additional issues appear; namely, the discrete equivalent Muskingum model should include the dynamics of a ZOH and second the Shannon criteria should be fulfilled when selecting the sampling time. These problems appear since, usually in hydraulics, the Muskingum model is discretized using the bilinear transformation and with a sampling time that provides hydrological meaning to the model (see Section3.1).However, when using such a model in closed-loop digital control, instead a ZOH transform should be used to take into account the dynamics of the ZOH that models the D/A conversion process. This is a standard procedure in classical digital control [54]. Then, equivalent discrete-time transfer function will have the same structure as that in the case of the bilinear transform, but the parameters of the model c
0
,
c
1, and c
2 will be different according to Table 1.If the discrete-time Muskingum model obtained by means of the bilinear transform is used to design the controller instead of the one obtained using the zero-order hold transform, the closed-loop behaviour will experiment a performance degradation when applied to a discrete-time Muskingum model using the zero-order hold transform according to the results presented in Figure7. The inclusion of the zero-order hold when applying digital control is unavoidable since the D/A converter is always present. As a conclusion, for closed-loop digital control design, the discrete-time Muskingum model should be obtained using the zero-order hold transform instead of the bilinear transform.Figure 7
Step closed-loop response for the case I: downstream flow with PI designed by continuous Muskingum with continuous Muskingum plant; case II: downstream flow with PI designed by discrete Muskingum (ZOH andT
s
=
60 s) with discrete Muskingum (ZOH and T
s
=
60 s) plant; case III: downstream flow with PI designed by discrete Muskingum (bilinear and T
s
=
600 s) with discrete Muskingum (bilinear and T
s
=
600 s) plant; case IV: downstream flow with PI designed by discrete Muskingum (bilinear and T
s
=
600 s) with discrete Muskingum (ZOH and T
s
=
600 s) plant; and case V: downstream flow with PI designed by discrete Muskingum (ZOH and T
s
=
600 s) with discrete Muskingum (ZOH and T
s
=
600 s) plant.


## 5.3. Effects of the Sampling Time

Even when discretizing the Muskingum model using the zero-order hold method, another problem appears when selecting the sampling time period satisfying the rule given in (9) in order that the discrete Muskingum model has hydrological meaning since the Shannon theorem [54] is not respected. Shannon’s theorem establishes that a good sampling time should take a value at least T
s
=
T
/
10, but according to rule (14), in the best case T
s
=
K
=
T
/
(
1
-
χ
). The selection of the sampling time following the Muskingum rule will derive again in a performance’s degradation according to Figure 7 with respect to the equivalent closed-loop continuous behaviour. In fact, in this figure it can be observed that the use of the sampling time suggested by Muskingum rule increases the overshoot. On the other hand, selecting the sampling time according to Shannon’s theorem derives in a closed-loop response that resembles that obtained when the continuous Muskingum model is used. However, even when using the ZOH transform method as a discretization procedure and the sampling time suggested by Shannon’s theorem, the problem of the structure in the Muskingum model described in Section 5.1 still remains.


## 6. Control Results for Real Plant

Finally, in order to validate both models (Muskingum and ID) for control design, the control results using the real plant (simulated by means of Saint-Venant equations) are presented in Figures8, 9, and 10. In this section, two usual methodologies for canal feedback control using a PI and with a Smith predictor (SP) scheme [55, 56] and predictive control [50, 57, 58] are tested for this study and for validation.Figure 8
Step closed-loop response for the case I: downstream flow with PI designed by discrete Muskingum (ZOH,T
s
=
60 s) with continuous Muskingum plant; case II: downstream flow with PI designed by discrete Muskingum (ZOH, T
s
=
600 s) with continuous Muskingum plant; case III: downstream flow with PI designed by discrete Muskingum (ZOH, T
s
=
60 s) with continuous ID plant; case IV: downstream flow with PI designed by discrete Muskingum (ZOH, T
s
=
600 s) with continuous ID plant; and case V: downstream flow with PI designed by continuous Muskingum with continuous Muskingum plant.Figure 9
Step closed-loop response of the real irrigation plant (Saint-Venant equations) using a PI1 and a PI2; step closed-loop response of the ID continuous plant using PI1 and step closed-loop response of the Muskingum continuous plant using PI2.Figure 10
Step closed-loop response of the real irrigation plant (Saint-Venant equations) using predictive control with a Muskingum model and an ID model (both with and without ZOH) for different predictive horizonλ and different T
s.


### 6.1. Using a PI Controller with a SP Scheme

For the PI with a SP scheme (to compensate the delay of the system) designed by pole placement technique [59–61], it can be observed that the results considering the real plant and the continuous ID model (used as a plant) using the PI designed with discrete ID model (ZOH, T
s
=
60 s) (denoted as PI1) are very similar. On the other hand, the closed-loop control results for the PI designed with discrete Muskingum model (bilinear, T
s
=
600 s) (denoted as PI2) considering the real plant are very different from the ones obtained using the continuous Muskingum model (used as a plant). These prove the validity of the ID model for control design, while showing that the Muskingum model is not suitable for this purpose.


### 6.2. Using a Predictive Controller with a SP Scheme

It is well known that predictive control essentially relies on the use of a model capable of forecasting the system output as a function of the system inputs on moving horizon setting and computing the control sequence making that predicted output could follow a desired trajectory through the minimization of a performance criterion. The predictive control is suitable for systems with delay and the model precision is not critical because a new predictive model is obtained from a new calibration every time the system operates.λ, the prediction horizon, is chosen in order to predict the transitory dynamics and to yield the adequate performance. The selection of prediction horizon will drive to a tradeoff between the abrupt control effort and the swift in the response. The precision and accuracy of predictive control have been studied with an ID model with T
s
=
60 s and Muskingum model with sample times that give physical meaning to Muskingum; that is, T
s
=
600 s. The ID model and Muskingum model provide a temporal response with approximately the same rapidity. However, the former presents more overshoot than the latter. Therefore, the predictive controller based on Muskingum model is more energetically efficient (less oscillations in the spillway movements). After analyzing the results obtained from the simulations, it is recommended to use the Muskingum model altogether with the predictive control strategy because it is a stationary, good prediction model and the delay is implicit.


## 6.1. Using a PI Controller with a SP Scheme

For the PI with a SP scheme (to compensate the delay of the system) designed by pole placement technique [59–61], it can be observed that the results considering the real plant and the continuous ID model (used as a plant) using the PI designed with discrete ID model (ZOH, T
s
=
60 s) (denoted as PI1) are very similar. On the other hand, the closed-loop control results for the PI designed with discrete Muskingum model (bilinear, T
s
=
600 s) (denoted as PI2) considering the real plant are very different from the ones obtained using the continuous Muskingum model (used as a plant). These prove the validity of the ID model for control design, while showing that the Muskingum model is not suitable for this purpose.


## 6.2. Using a Predictive Controller with a SP Scheme

It is well known that predictive control essentially relies on the use of a model capable of forecasting the system output as a function of the system inputs on moving horizon setting and computing the control sequence making that predicted output could follow a desired trajectory through the minimization of a performance criterion. The predictive control is suitable for systems with delay and the model precision is not critical because a new predictive model is obtained from a new calibration every time the system operates.λ, the prediction horizon, is chosen in order to predict the transitory dynamics and to yield the adequate performance. The selection of prediction horizon will drive to a tradeoff between the abrupt control effort and the swift in the response. The precision and accuracy of predictive control have been studied with an ID model with T
s
=
60 s and Muskingum model with sample times that give physical meaning to Muskingum; that is, T
s
=
600 s. The ID model and Muskingum model provide a temporal response with approximately the same rapidity. However, the former presents more overshoot than the latter. Therefore, the predictive controller based on Muskingum model is more energetically efficient (less oscillations in the spillway movements). After analyzing the results obtained from the simulations, it is recommended to use the Muskingum model altogether with the predictive control strategy because it is a stationary, good prediction model and the delay is implicit.


## 7. Conclusions

In this paper, the use of the Muskingum model in closed-loop digital and continuous control has been analyzed and discussed regarding three main issues: the structure of the model, the discretization method, and the selection of the sampling time. Regarding the structure of the Muskingum model in continuous time, it can be observed that it introduces a nonminimum-phase zero that produces negative responses at the beginning of the response to a step input. This behaviour is not present in the real canal making this model not useful in continuous time. Such nonminimum phase behaviour is avoided in discretizing the Muskingum model using the rule that guarantees the sampling time with the order of the delay. However, such a selection of the sampling time is not supported by Shannon’s theorem. Moreover, the discretization method typically used with the Muskingum model is the bilinear transform that does not take into account the zero-order hold. However, even when using the zero-order hold method as a discretization procedure and sampling time fulfilling Shannon’s theorem, the problem of the nonminimum phase introduced by the unstable zero appears again since the closed-loop behaviour in this case resembles the equivalent closed-loop continuous behaviour. As a conclusion, the Muskingum model is a good open-loop prediction model in discrete time using the suggested discretization rule in order to avoid the nonminimum phase behaviour and it is suitable for flood routing applications. For control purposes, mathematically simple control laws that can be designed (either based on the PID law or based on a predictive control strategy) with both models have been considered. These laws should allow an easy and swift computation, guaranteeing stability and convergence to the desired setpoint. However, when using the Muskingum model for designing a closed-loop controller, the results are not as good as expected. In particular, if the controller is implemented in a digital way, the mandatory inclusion of the ZOH effect and selection of the sampling time according to Shannon’s theorem derive in the reappearance of the nonminimum phase behaviour that Muskingum sampling time rule avoided. Furthermore, if the sampling time is reduced in the simulation of Muskingum controller, there is an increment of oscillations in the closed-loop response. This makes the Muskingum model also not useful for digital feedback control design. But in the case of predictive control, the explained Muskingum behaviour does not influence the predictive control algorithm because this type of control requires only a good prediction model at each prediction horizon. Since the delay is implicit in the system, a prediction model is obtained when a suitable sample time is selected. Therefore, with the use of Muskingum model, energetically efficient control results are obtained, although ID model also presents reasonable results. On the other hand, when PID controllers are designed using the ID model, the control results are temporally more efficient (less time of convergence and more rapidity) and more energetically efficient than Muskingum model. However, Muskingum model is more energetically efficient when considering a predictive law, although the speed of the response of control results is similar using both models.



---

*Source: 197907-2014-08-14.xml*