Urban Comprehensive Water Consumption: Nonlinear Control of Production Factor Input Based upon the C-D Function

Abstract

Utilizing the urban water demand function and the Cobb-Douglas (C-D) production function, an economic control model for the multi-input-multi-output (MIMO) nonlinear system was designed and implemented to describe urban comprehensive water consumption, where the urban water demand function was expressed as the product of the number of water users and per capita comprehensive water consumption, and the urban water supply function was expressed as a C-D production function. The control variables included capital investment and labor input for the urban water supply. In contrast to the Solow model, Shell model and aggregate model with renewable labor resources, the proposed model eliminated value constraints on investment and labor input in the state equations and hence avoided the difficulty in applying these models to urban water supply institutions. Furthermore, the feedback linearization control design (FLCD) method was employed to accomplish stability of the system. In contrast to the optimal control method, the FLCD method possesses an explicit solution of the control law and does not require the solution of a two-point boundary value problem of an ordinary differential equation, making the method more convenient for application. Moreover, two different scenarios of urban water consumption, one for the growth period and the other for the decline period, were simulated to demonstrate the effectiveness of the proposed control scheme.

1. Introduction

China has a great shortage of water resources. According to the world development indicators issued by the World Bank in 2017, China’s renewable internal freshwater resources per capita was 2062 m³ in 2014, which was only 35% of the world average of 5917 m³ in that same year [1]. In 2000, the Chinese Academy of Engineering predicted that China’s per capita water resources will have fallen to 1760m³ by 2030 [2].

In order to alleviate the shortage of water resources, China’s Central Government issued the “Decision on Accelerating the Reform and Development of Water Conservancy” in 2011, administratively ordering an increase of water conservancy investment nationwide to a level of 4 trillion Yuan during 2011–2020. This means that the average annual investment is 400 billion Yuan, which is twice as high as it was in 2010.

Consequently, the input of public finance at all levels of the water conservancy industry has increased significantly in both the total amount and the increase rate. However, how to adjust this input, which is up to 4 trillion Yuan, to water demand over a period of time is a question worth studying. The remainder of the paper is dedicated to this issue.

As a widely used production function, the Cobb-Douglas (C-D) production function also has some applications in resource management and water management.

In 2012, Kiselev et al. considered a resource distribution problem on a finite time interval with a terminal function for a two-sector economic model with a two-factor C-D production function with distinct amortization factors [3]. In 2013, Zhang et al. established a regional water demand prediction model based on the C-D model, through which the contribution rates of regional water demand influencing factors could be calculated [4]. In 2015, Li used the product of logistic equation and per capita water demand to denote domestic water demand, and the C-D production function to denote water supply. Through the dynamic coordination of capital and labor input, the urban water supply dynamic was made to meet domestic water demand [5]. In 2017, Zhang et al. developed a fuzzy credibility-constrained interval two-stage stochastic programming (FCITSP) model based on the C-D function to optimize water distribution of three industries based on a water demand prediction under multiple uncertainties [6].

Furthermore, there have been some impressive achievements in the application of the dynamical control theory in water management.

In 2001, Mourad et al. studied the discrete-time adaptive control method of water supply systems by controlling the flow rate of water pumps to regulate water consumption [7]. In 2003, Eker et al. adopted the robust polynomial $H_{\infty }$ method to achieve the control of a water supply simulation system consisting of pumping stations and intermediate storage reservoirs [8]. In 2007, Izquierdo et al. combined deterministic and machine learning models to establish a complex hybrid model for water system fault detection [9]. In 2010, 2011 and 2015, Li et al. studied urban water supply control using the optimal control method [5,10,11], and in 2010 and 2018, Li et al. studied robust control of urban water supply based on the pole allocation method [12,13] and the Lyapunov stability theory combined with the linear matrix inequality method [14]. In 2015, Sun et al. studied the model predictive control method for the optimal allocation of large-scale real-time water supply systems and applied it to the water supply network in Catalonia, Spain [15]. In 2015, Sarbu et al. optimized the water supply of single-source or multi-source circular distribution networks according to change of demand [16]. Finally, in 2016, Pereira et al. studied the application of a robust predictive controller to track drinking water networks with periodicity in Barcelona, ​​Spain [17].

The focus of the existing literature has been on physical control of water pumps or tanks [7,8,9,15,16,17]; the optimal control method [5,10,11] is a commonly used technique in dynamic system control, which can get the system dynamic optimal solution in a certain period of time. However, this method has some difficulties in its applications: It needs to solve a two-point boundary value problem of differential equations involving state variables and covariate variables [18,19].

Moreover, robust control based on the pole placement method [12,13] needs to build a complex system compensator, which is difficult to implement. Additionally, the robust control method [14] only considers linear systems.

The existing literature has paid little attention to the following situation: The nonlinear system of urban comprehensive water consumption is composed of the water demand function and the water supply function, where the former is described based on per capita comprehensive water consumption, and the latter is described by the C-D production function. This comprises the main study of this paper.

We elaborate on and outline the study below.

The goal of system management is to ensure the stability of supply and demand of urban comprehensive water consumption. In order for us to achieve this primary goal on the system, a nonlinear economic control method based on feedback linearization was employed in our model. Through dynamic control of production factors, i.e., capital investment and labor input, water supply can meet dynamic water demand.

At present, this kind of research is rare in the literature; however, it is of great significance to the optimization of input efficiency of the urban water industry.

From an economic management point of view, this paper studies production factor input control of urban comprehensive water consumption systems. Based on per capita comprehensive water consumption and the C-D production function, a multi-input-multi-output (MIMO) nonlinear control system model of urban comprehensive water consumption was developed. We used capital investment and labor input of urban water supply as control variables. The model developed in this paper was compared with some economic dynamic models that are used in economic, ecological and environmental sciences, such as the Solow model, the Shell model and the aggregate model with renewable labor resources. The state equations of capital stock and labor stock of this proposed model were different from the aforementioned models as this new model canceled value constraints of these models on investment and labor input and avoided the situation where these models are difficult to apply to urban water supply agencies.

In order to avoid difficulties in the application of the optimal control method, this paper adopted the feedback linearization control design method to achieve the stable supply and demand of the urban comprehensive water consumption system. In contrast to the optimal control method, the feedback linearization control method can get an explicit solution of the control law, but it does not need to solve the complex two-point boundary value problems of differential equations, which makes the feedback linearization control method more convenient for solution and application. The simulation results showed the effectiveness of the control schemes in different situations.

Through feedback linearization control of Equation (22), the urban comprehensive water consumption nonlinear control system (Equation (6)) can be transformed into the linear feedback control system (Equation (19)).

In the remainder of the paper, Section 2 introduces the economic control model for urban comprehensive water consumption. Section 3 compares the proposed model with the existing economic dynamic models. In Section 4, we present the control design for the input of production factors regarding the urban comprehensive water consumption system. Section 5 illustrates the proposed model with a case simulation. For the theoretical solution, we further simulated the control system with MATLAB by assuming two different situations on urban comprehensive water consumption: Increasing period and decreasing period with respect to different feedback gain values. Section 6 discovers the explicit solution based on the feedback linearization control method. Furthermore, a simulation with Excel is performed in this section. Finally, in Section 7, the conclusion is drawn.

2. Economic Control Model for Urban Comprehensive Water Consumption

Let $D\left(t\right)$ be the demand function for urban comprehensive water consumption, where $t$ stands for time. This can be expressed as the product of the number of urban water users, $M\left(t\right)$, and per capita comprehensive water consumption, $b$, i.e.,

(1)$D\left(t\right)=bM\left(t\right).$

This is a widely used urban water demand prediction method [20].

Let $S\left(t\right)$ be the supply function of urban comprehensive water consumption. From the perspective of economic management, this is mainly a function of capital and labor. Therefore, according to the theory of water resource economics [21], it can be expressed by the C-D production function, which is commonly used in economics. The C-D function has been widely used to measure the relationship between input and output, including in the field of water management [4,5,21]. We have

(2)$S\left(t\right)=AK{\left(t\right)}^{\alpha }L{\left(t\right)}^{\beta },$

where $K\left(t\right)$ is the capital stock of urban water supply agencies, $L\left(t\right)$ is the labor stock of urban water supply agencies, $A$ is the technical level coefficient of urban water supply agencies, $A>0$, $\alpha$ is the elasticity coefficient of the capital output, and $0<\alpha <1$, and $\beta$ is the elasticity coefficient of the labor output, $0<\beta <1$.

The state equation of $K\left(t\right)$ can be expressed as

(3)$\dot{K}\left(t\right)=I\left(t\right)-\delta K\left(t\right),$

where $I\left(t\right)$ is the amount of capital investment by urban water supply agencies and $\delta$ is the depreciation rate of capital, $0<\delta <1$.

The state equation of $L\left(t\right)$ can be expressed as

(4)$\dot{L}\left(t\right)=P\left(t\right)-\sigma L\left(t\right),$

where $P\left(t\right)$ is the labor input of urban water supply agencies and $\sigma$ is the natural withdrawal rate of labor force due to retirement or voluntary resignations, $0\le \sigma <1$.

It is assumed that the available labor in urban water supply agencies is proportional to the population of urban water users, namely

(5)$L\left(t\right)=aM\left(t\right),$

where $a$ is the ratio coefficient of the stock of labor in urban water supply agencies to the population of urban water users, $0<a<1$.

If we rearrange the above equations, by substituting Equation (5) with Equation (1), and combining this with Equations (2)–(4) we have the following model:

(6)$\{\begin{array}{l}D\left(t\right)=\frac{b}{a}L\left(t\right),0<a<1\\ S\left(t\right)=AK{\left(t\right)}^{\alpha }L{\left(t\right)}^{\beta },0<\alpha <1,0<\beta <1\\ \dot{K}\left(t\right)=I\left(t\right)-\delta K\left(t\right),0<\delta <1\\ \dot{L}\left(t\right)=P\left(t\right)-\sigma L\left(t\right),0\le \sigma <1\end{array}$

In the above model, the last two differential equations are linear in $K\left(t\right)$ and $L\left(t\right)$, but in the second algebraic equation, the C-D production function $S\left(t\right)$ is nonlinear in $K\left(t\right)$ and $L\left(t\right)$.

Since $S\left(t\right)$ is nonlinear in $K\left(t\right)$ and $L\left(t\right)$, the urban comprehensive water consumption system described in Equation (6) is a MIMO nonlinear control system with two dimensional inputs and two dimensional outputs, where $I\left(t\right)$ and $P\left(t\right)$ are two-dimensional control input variables, $D\left(t\right)$ and $S\left(t\right)$ are two-dimensional output variables, and $K\left(t\right)$ and $L\left(t\right)$ are two-dimensional state variables.

The control objectives of the above model are as follows: The system output $D\left(t\right)$ and $S\left(t\right)$ will, over time, tend to and maintain the desired outputs $D*$ and $S*$; that is, water demand and supply will tend to each ideal value respectively, and system stability will be realized through the following optimization problem (Equation (8)), where the objective functions are given by

(7)$\{\begin{array}{l}{J}_1=\min \left\{{\left[\underset{t}{\lim }D\left(t\right)-D*\right]}^2\right\}\\ J_2=\min \left\{{\left[\underset{t}{\lim }S\left(t\right)-S*\right]}^2\right\}\end{array}.$

In general, $S*=D*$; that is, at a certain point in time the system realizes the balance of supply and demand. According to individual needs one can assume $S*\ne D*$, such as a setting of $S*>D*$, to keep urban water supply capacity to a certain degree of redundancy.

By combining this with Equations (6) and (7), we have the nonlinear control model of production factor input in a urban comprehensive water consumption model based on the C-D function:

(8)$\{\begin{array}{l}{J}_1=\min \left\{{\left[\underset{t}{\lim }D\left(t\right)-D*\right]}^2\right\}\\ J_2=\min \left\{{\left[\underset{t}{\lim }S\left(t\right)-S*\right]}^2\right\}\\ s.t.\\ \{\begin{array}{l}D\left(t\right)=\frac{b}{a}L\left(t\right),0<a<1\\ S\left(t\right)=AK{\left(t\right)}^{\alpha }L{\left(t\right)}^{\beta },0<\alpha <1,0<\beta <1\\ \dot{K}\left(t\right)=I\left(t\right)-\delta K\left(t\right),0<\delta <1\\ \dot{L}\left(t\right)=P\left(t\right)-\sigma L\left(t\right),0\le \sigma <1\end{array}\end{array}$

The goal is to choose $I\left(t\right)$ and $P\left(t\right)$ to simultaneously minimize $J_1$ and $J_2$.

Using the above model, from the perspective of economic management, we hope to achieve the goal of stabilizing urban comprehensive water consumption systems through economic control of the two production factors, $I\left(t\right)$ and $P\left(t\right)$.

3. Consumption Compared with the Existing Economic Dynamics Models

Currently, there appears to be no dynamic model for urban comprehensive water consumption with the same control variables and state equations as described in Section 2. Our model can be compared with some similar models used in the economic, ecological and environmental sciences, such as the Solow model, the Shell model and the aggregate model with renewable labor resources [22].

For example, the Solow model is

(9)$\{\begin{array}{l}I\left(t\right)=qS\left(t\right),0<q<1\\ \dot{K}\left(t\right)=I\left(t\right)-\delta K\left(t\right),0<\delta <1\\ \dot{L}\left(t\right)=\eta L\left(t\right)\end{array}.$

In this case, $q$ is the control variable of the model by fixing the ratio of $I\left(t\right)$ and $S\left(t\right)$ to a constant value. Due to $0<q<1$, this also means that $I\left(t\right)$ cannot take a negative value, and capital investment cannot be withdrawn. Therefore, the value of $I\left(t\right)$ is strongly constrained. In addition, $\eta$, in the state equation of $L\left(t\right)$, is the relative rate of change of the labor force, which is a constant and results in exponential growth. Finally, $\dot{L}\left(t\right)$ does not have the control variable $P\left(t\right)$.

The Solow Model (Equation (9)) becomes the Shell Model (Equation (10)) when q is allowed to become a function of time. The Shell model becomes:

(10)$\{\begin{array}{l}I\left(t\right)=q\left(t\right)S\left(t\right),0\le q\left(t\right)\le 1\\ \dot{K}\left(t\right)=I\left(t\right)-\delta K\left(t\right),0<\delta <1\\ \dot{L}\left(t\right)=\eta L\left(t\right)\end{array}$

Although this relaxes some constraints on $I\left(t\right)$, the constraint that $I\left(t\right)$ cannot take a negative value still exists. Furthermore, there is still no control variable $P\left(t\right)$ in $\dot{L}\left(t\right)$.

In order to apply $P\left(t\right)$ to control $L\left(t\right)$, another modification of the Solow Model (Equation (9)) is the aggregate model with renewable labor resources:

(11)$\{\begin{array}{l}I\left(t\right)=q_{1}S\left(t\right),0<q_1<1\\ \dot{K}\left(t\right)=I\left(t\right)-\delta K\left(t\right),0<\delta <1\\ P\left(t\right)=q_{2}S\left(t\right),0<q_2<1\\ \dot{L}\left(t\right)=P\left(t\right)-\eta L\left(t\right)\end{array}$

Compared with the Solow Model (Equation (9)), the aggregate model (Equation (11)) adds a new parameter to the model: $P\left(t\right)$, which is used to reproduce the labor force. The model sets $P\left(t\right)$ to a constant fixed ratio $q_2$ of $S\left(t\right)$. Therefore, $q_1$ and $q_2$ become two independent and unknown control variables. Similar to $q_1$ in $I\left(t\right)$, $q_2$ satisfies $0<q_2<1$, which means that $P\left(t\right)$ cannot take a negative value in the aggregate model (Equation (11)).

The constraints, listed above, on $I\left(t\right)$ or $P\left(t\right)$ in the Solow Model and the aggregate model are difficult to apply by urban water agencies. Because urban water supply agencies may pursue the goal of a balance between supply and demand, a situation in which water supply is greater than water demand means that $I\left(t\right)$ or $P\left(t\right)$ may need to take a negative value.

The model established in this paper (Equation (8)), eliminates the constraints on $I\left(t\right)$ and $P\left(t\right)$. This means that urban water supply agencies can determine the appropriate investment and labor input at each stage of $t$ in accordance with the actual situation, so that the system goals can be achieved better.

4. Control Design of the Input of Production Factors for Urban Comprehensive Water Consumption Systems

In the analysis and design of a linear system, for the convenience of notations and analysis, the system is often transformed so that its equilibrium point becomes the zero point of the state space. Such a transformation can also be made for a specific equilibrium point for the non-linear system:

(12)$\dot{x}=f\left(x\right).$

If the equilibrium points to be studied is $x*$, and a new variable $y$ is introduced,

$y=x-x*$

the equation $y$ is obtained by substituting $x=y+x*$ in Equation (12):

(13)$\dot{y}=f\left(y+x*\right)$

It is easy to verify that the solution of Equation (13) corresponds to that of Equation (12), and that $y=0$ corresponds to $x=x*$, which is an equilibrium point of Equation (13). Therefore, if we want to study the behavior of Equation (12) near the equilibrium point $x*$, we study precisely the behavior of Equation (13) in the neighborhood of the origin [23].

Therefore, in our model (Equation (8)), we first investigate the expected output $D*=0$ and $S*=0$; that is, the equilibrium point of the system is the origin. Then, the above transformation is used to investigate the situation of $D*$ and $S*$ for any other point.

For the case of the equilibrium point as the original point, the feedback linearization method is used to control the design in the following way.

Again, the goal of the system control process is to achieve stability for the urban comprehensive water consumption system by selecting appropriate values for the two control variables $I\left(t\right)$ and $P\left(t\right)$.

First, it is necessary to determine whether the system has internal dynamics, which can be accomplished by solving the total relative degree of the system. For this reason, the output variables $D\left(t\right)$ and $S\left(t\right)$ in Equation (8) are differentiated until the input $I\left(t\right)$ or $P\left(t\right)$ appears. We have

(14)$\{\begin{array}{l}\dot{D}=\frac{b}{a}\dot{L}=\frac{b}{a}\left(P-\sigma L\right)\\ \dot{S}=A\alpha K^{\alpha -1}{L}^{\beta }\dot{K}+A\beta K^{\alpha }{L}^{\beta -1}\dot{L}=A\alpha K^{\alpha -1}{L}^{\beta }\left(I-\delta K\right)+A\beta K^{\alpha }{L}^{\beta -1}\left(P-\sigma L\right)\end{array}$

From Equation (14), we can see that the input $I\left(t\right)$ or $P\left(t\right)$ has appeared in $\dot{D}\left(t\right)$ and $\dot{S}\left(t\right)$ when $D\left(t\right)$ and $S\left(t\right)$ are differentiated. Therefore, Equation (8) has the relative degree of $\left(1,1\right)$ at the origin, and the total relative degree of the system is $r=1+1=2$. The total relative degree is equal to the dimension of the system state variable, and hence the system has no internal dynamics.

Equation (14) is rewritten as a matrix form,

(15)$\left(\begin{array}{c}\dot{D}\\ \dot{S}\end{array}\right)=\left(\begin{array}{c}-b\sigma L/a\\ -\left(\alpha \delta +\beta \sigma \right)A{K}^{\alpha }{L}^{\beta }\end{array}\right)+\left(\begin{array}{cc}0& b/a\\ A\alpha K^{\alpha -1}{L}^{\beta }& A\beta K^{\alpha }{L}^{\beta -1}\end{array}\right)\left(\begin{array}{c}I\\ P\end{array}\right)$

Let

$E=\left(\begin{array}{cc}0& b/a\\ A\alpha K^{\alpha -1}{L}^{\beta }& A\beta K^{\alpha }{L}^{\beta -1}\end{array}\right)$

and we can see that $E$ is an invertible matrix. If the control input transformation is

(16)$u=\left(\begin{array}{c}I\\ P\end{array}\right)=E^{-1}\left(\begin{array}{c}{v}_1+b\sigma L/a\\ v_2+\left(\alpha \delta +\beta \sigma \right)A{K}^{\alpha }{L}^{\beta }\end{array}\right),$

$v_1,v_2$ are the new undetermined inputs. From Equations (15) and (16), we have

(17)$\{\begin{array}{l}\dot{D}=v_1\\ \dot{S}=v_2\end{array}.$

From Equation (17), $v_1,v_2$ are only applied to $D\left(t\right)$ and $S\left(t\right)$, respectively, so Equation (16) is the decoupling control law, and the invertible matrix $E$ is the decoupling matrix of the system. Furthermore, because the system has no internal dynamics, the input-state linearization of the nonlinear system (Equation (8)) can be realized by using the control law of Equation (16). Therefore, using the equivalent input of Equation (17), which is similar to that of a single-input-single-output system, we can solve the system’s stabilization and tracking problems without having to consider the stability of internal dynamics.

Let the new input be

(18)$\{\begin{array}{l}{v}_1=-k_{1}D\\ v_2=-k_{2}S\end{array}$

where, $k_1$ and $k_2$ are undetermined feedback gains. By substituting the state feedback of Equation (18) with Equation (17), the system in Equation (17) is then transformed into the following linear feedback control system:

(19)$\{\begin{array}{l}\dot{D}=-k_{1}D\\ \dot{S}=-k_{2}S\end{array}$

For the system in Equation (19), by just selecting $k_1>0,k_2>0$, one can guarantee the exponential stability of $D\left(t\right)$ and $S\left(t\right)$ and realize the equilibrium of the system at the origin. The solution of the system in Equation (19) can be given in an explicit form:

(20)$\{\begin{array}{l}D=D\left(0\right)e^{-k_{1}t}\\ S=S\left(0\right)e^{-k_{2}t}\end{array}.$

By substituting Equation (18) with Equation (16), the value of the original control input variable is obtained as

(21)$u=\left(\begin{array}{c}I\\ P\end{array}\right)=E^{-1}\left(\begin{array}{c}-k_{1}D+b\sigma L/a\\ -k_{2}S+\left(\alpha \delta +\beta \sigma \right)A{K}^{\alpha }{L}^{\beta }\end{array}\right)=E^{-1}\left(\begin{array}{c}-k_{1}D+b\sigma L/a\\ \left(-k_2+\alpha \delta +\beta \sigma \right)S\end{array}\right)$

Equation (21) completes the control design of Equation (8) at the origin.

The following is the control design when the equilibrium points $D*$ and $S*$ of the output of Equation (8) are at any other point.

In this case, since Equations (12) and (13) have a one-to-one correspondence relationship, we can take the value of the control input variable of Equation (21) as

(22)$u=\left(\begin{array}{c}I\\ P\end{array}\right)=E^{-1}\left(\begin{array}{c}-k_1\left(D-D*\right)+b\sigma L/a\\ \left(-k_2+\alpha \delta +\beta \sigma \right)\left(S-S*\right)\end{array}\right)$

Equation (22) gives the values of $I\left(t\right)$ and $P\left(t\right)$. Accordingly, we can determine the input of the production factor for the urban comprehensive water consumption system, and then the control goal can be realized when the equilibrium points $D*$ and $S*$ of the system output are at any other point.

For the model established in this paper, it can be seen from the simulation below that the key to determine Equation (22) is to select the appropriate $k_1$ and $k_2$, which have an important impact on the speed of the realization of the system goal.

5. Case Simulation

In accordance with the economic evaluation method and parameter of a construction project [24], the parameters of the urban comprehensive water consumption system were set as follows: $b=80$ m³, $a=0.002$, $\alpha =0.7$, $\beta =0.3$, $A=20$, $\delta =0.04$, $\sigma =0.02$, $t_0=0$, $t_f=20$ a, where $t_f$ was the simulation termination time and a was the number of years, $K\left(0\right)=1500$ × 10⁴ Yuan, and $L\left(0\right)=0.1$ × 10⁴ people. From these we got $D\left(0\right)=4000$ × 10⁴ m³ and $S\left(0\right)=1676.1$ × 10⁴ m³. The control goal was 20 years. By adjusting $I(t)$ and $P(t)$, the system outputs $D\left(t\right)$ and $S\left(t\right)$ could reach and maintain the expected outputs $D*$ and $S*$, respectively. The condition, $S*=D*$, could achieve the balance between supply and demand.

The following two cases of $D*$ were investigated.

Case 1: $D*=5000$ × 10⁴ m³. There was $D*>D(0)$ at this time, indicating that the urban comprehensive water consumption system was in a period of rising.

Case 2: $D*=3000$ × 10⁴ m³. There was $D*<D(0)$ at this time, indicating that the urban comprehensive water consumption system was in a period of decline.

The World Water Vision made a prediction on urban water consumption, stating that by 2025, it will have risen rapidly in developing countries and will have stayed stable or declined in developed countries [25]. Therefore, the simulation of the above two cases of $D*$ could correspond to the two predictions of the World Water Vision.

Both of these cases were simulated in the following way.

For the feedback gains $k_1$ and $k_2$ in Equation (18), the conditions of $D\left(t\right)$ and $S\left(t\right)$ exponential stability can be satisfied only by $k_1>0,k_2>0$. However, different values of $k_1$ and $k_2$ have a different influence on system states. For the values of $k_1$ and $k_2$, it is necessary, through simulation, to determine the speed requirement of achieving the system target and the actual value range of the control variables.

Following the simulation of some value combinations of $k_1$ and $k_2$, the effects of different combinations of the values on the above two cases of the system were investigated. Taking $k_1$ as 0.1, 0.9 and 1.8 and $k_2$ as 1, 1.8 and 2.6, their combinations were simulated, respectively. There were a total of 9 combinations:

(1) $k_1=0.1,k_2=1$, (2) $k_1=0.1,k_2=1.8$, (3) $k_1=0.1,k_2=2.6$, (4) $k_1=0.9,k_2=1$, (5) $k_1=0.9,k_2=1.8$, (6) $k_1=0.9,k_2=2.6$, (7) $k_1=1.8,k_2=1$, (8) $k_1=1.8,k_2=1.8$, (9) $k_1=1.8,k_2=2.6$.

After substituting Equation (22) with Equation (15), the system was transformed into a linear differential equation by feedback linearization control. The initial value problem of the first-order differential equation was

(23)$\{\begin{array}{l}\dot{x}=f\left(t,x\right),t\in \left[t_0,t_f\right]\\ x\left(t_0\right)=x_0\end{array}.$

The fourth-order and fifth-order Runge-Kutta-Fehlberg algorithm with a variable step size can be used to solve iteratively [26].

• ①. Set the value of $t_0,t_f,x_0,h=0.1,\epsilon =1\mathrm{e}-7$, where $h$ is the initial step size, and $\epsilon$ is the specified error control tolerance. Let $n=0$.

• ②. Let $t_{n+1}=t_n+h$. If $t_{n+1}>t_f$, go to ⑤; otherwise calculate

  $\{\begin{array}{l}{\tilde{x}}_{n+1}=x_n+h\left(\frac{25}{216}{K}_1+\frac{1408}{2565}{K}_3+\frac{2197}{4104}{K}_4-\frac{1}{5}{K}_5\right)\\ {\widehat{x}}_{n+1}=x_n+h\left(\frac{16}{135}{K}_1+\frac{6656}{12825}{K}_3+\frac{28561}{56430}{K}_4-\frac{9}{50}{K}_5+\frac{2}{55}{K}_6\right)\end{array}$

  where

  $\{\begin{array}{l}{K}_1=f\left(t_n,x_n\right)\\ K_2=f\left(t_n+\frac{h}{4},x_n+\frac{h}{4}{K}_1\right)\\ K_3=f\left(t_n+\frac{3h}{8},x_n+\frac{3h}{32}{K}_1+\frac{9h}{32}{K}_2\right)\\ K_4=f\left(t_n+\frac{12h}{13},x_n+\frac{1932h}{2197}{K}_1-\frac{7200h}{2197}{K}_2+\frac{7296h}{2197}{K}_3\right)\\ K_5=f\left(t_n+h,x_n+\frac{439h}{216}{K}_1-8h{K}_2+\frac{3680h}{513}{K}_3-\frac{845h}{4104}{K}_4\right)\\ K_6=f\left(t_n+\frac{h}{2},x_n-\frac{8h}{27}{K}_1+2h{K}_2-\frac{3544h}{2565}{K}_3+\frac{1859h}{4104}{K}_4-\frac{11h}{40}{K}_5\right)\end{array}$

• ③. Calculate

  $\frac{\left|{\widehat{x}}_{n+1}-{\tilde{x}}_{n+1}\right|}{h}=\left|\frac{1}{360}{K}_1-\frac{128}{4275}{K}_3-\frac{2197}{75240}{K}_4+\frac{1}{50}{K}_5+\frac{2}{55}{K}_6\right|$

  If $\left|{\widehat{x}}_{n+1}-{\tilde{x}}_{n+1}\right|/h<\epsilon$, then output $t_{n+1},{\tilde{x}}_{n+1}$, let $n=n+1$, and go to ②; otherwise go to ④.

• ④. If

  $\{\begin{array}{l}q={\left(\frac{\epsilon h}{2\left|{\widehat{x}}_{n+1}-{\tilde{x}}_{n+1}\right|}\right)}^{1/4}\\ h=qh\end{array}$

  go to ②.

• ⑤. End.

The fourth-order and fifth-order Runge-Kutta-Fehlberg algorithm with a variable step size can be implemented by the differential equation solving function ode45( ) of MATLAB software [27], which can be used in this simulation. To improve the calculation accuracy of the simulation results, the default value of the relative error tolerance RelTol in the control parameter was changed from 1e-3 to 1e-7.

The simulation results of urban comprehensive water consumption system control under the above two cases were calculated, as shown in Figure 1, Figure 2, Figure 3, Figure 4, Figure 5, Figure 6, Figure 7, Figure 8 and Figure 9. The left part of each figure corresponds to the system state of Case 1 and the right part corresponds to the system state of Case 2.

The simulation results show that as $k_1$ and $k_2$ increased, the range of $P(t)$ and $I(t)$ increased, and the control efforts in the short term also increased. The increase in control efforts meant that the target of the system could be realized faster, but demand for the input of the production factor was higher in the short term. Therefore, the specific value of $k_1$ and $k_2$ should be determined depending on the actual situation.

It is necessary to consider both the rising and the decline cases of urban comprehensive water consumption systems in a comprehensive way. In Combinations (1)–(3), because the value of $k_1$ was small, in Figure 1, Figure 2 and Figure 3, it was found that $S(t)>D(t)$ for Case 1 for a long time, and for Case 2, it was $S(t)<D(t)$. In Combination (4), although the value of $k_1$ increased, the value of $k_2$ was small, which caused, in Figure 4, $S(t)<D(t)$ in Case 1. The cases of Combinations (1)–(4) are not expected in reality, and these combinations can be discarded. In Combinations (5)–(9), because the values of $k_1$ and $k_2$ were larger, the input of $P(t)$ or $I(t)$ increased, $S(t)$ could approach $D*$ and achieve system stability in Figure 5, Figure 6, Figure 7, Figure 8 and Figure 9, and hence all these combinations could be selected according to actual conditions.

The simulation results in Case 1 and Case 2 show that the above feedback linearization control design scheme can adapt to the needs of the control of production factor input for urban comprehensive water consumption systems under varying conditions, coordinate water demand and water supply to an ideal value, and achieve the goal of system stability.

6. Explicit Solution Based on Feedback Linearization Control

The focus of this paper was to determine the proper value of the control variable, i.e., Equation (22). Through this formula, the nonlinear system can be converted into a linear system for control. The case simulation in the previous section was calculated directly using Equation (22).

Furthermore, it is also possible to find the explicit solution for the system that corresponds to Equation (22), and hence Excel can be used for a direct calculation. The following is the calculation process.

Analogous to the correspondence of Equation (21) to Equation (20), Equation (22) corresponds to Equation (24) below:

(24)$\{\begin{array}{l}\left(D-D*\right)=\left(D\left(0\right)-D*\right)e^{-k_{1}t}\\ \left(S-S*\right)=\left(S\left(0\right)-S*\right)e^{-k_{2}t}\end{array}$

i.e.,

(25)$\{\begin{array}{l}D=\left(D\left(0\right)-D*\right)e^{-k_{1}t}+D*\\ S=\left(S\left(0\right)-S*\right)e^{-k_{2}t}+S*\end{array}$

By Equation (25), we have

(26)$\{\begin{array}{l}\dot{D}=-k_1\left(D\left(0\right)-D*\right)e^{-k_{1}t}\\ \dot{S}=-k_2\left(S\left(0\right)-S*\right)e^{-k_{2}t}\end{array}$

From Equations (1) and (5), it holds that

(27)$L=\frac{a}{b}D$

By Equations (1), (4) and (5), we obtain

(28)$\dot{D}=\frac{b}{a}\left(P-\sigma L\right)$

From Equations (25)–(28), we have

(29)$P=\left(\sigma -k_1\right)\left(\frac{a}{b}\right)\left(D\left(0\right)-D*\right)e^{-k_{1}t}+\sigma \left(\frac{a}{b}\right)D*$

By Equation (2), we get

(30)$K={\left(\frac{1}{A}S{L}^{-\beta }\right)}^{\frac{1}{\alpha }}$

and

(31)$\dot{S}=S\left(\alpha \frac{\dot{K}}{K}+\beta \frac{\dot{L}}{L}\right)$

Equation (31) gives

(32)$\dot{K}=\frac{K}{\alpha }\left(\frac{\dot{S}}{S}-\beta \frac{\dot{L}}{L}\right)$

From Equations (32) and (3), we have

(33)$I=\frac{K}{\alpha }\left(\frac{\dot{S}}{S}-\beta \frac{\dot{L}}{L}\right)+\delta K.$

Combining Equations (25)–(27), (29), (30) and (33), we have obtained the following:

(34)$\{\begin{array}{l}D=\left(D\left(0\right)-D*\right)e^{-k_{1}t}+D*\\ S=\left(S\left(0\right)-S*\right)e^{-k_{2}t}+S*\\ \dot{D}=-k_1\left(D\left(0\right)-D*\right)e^{-k_{1}t}\\ \dot{S}=-k_2\left(S\left(0\right)-S*\right)e^{-k_{2}t}\\ L=\frac{a}{b}D\\ K={\left(\frac{1}{A}S{L}^{-\beta }\right)}^{\frac{1}{\alpha }}\\ I=\frac{K}{\alpha }\left(\frac{\dot{S}}{S}-\beta \frac{\dot{L}}{L}\right)+\delta K\\ P=\left(\sigma -k_1\right)\left(\frac{a}{b}\right)\left(D\left(0\right)-D*\right)e^{-k_{1}t}+\sigma \left(\frac{a}{b}\right)D*\end{array}$

Equation (34) gives the explicit solution for the system that corresponds to the control law of Equation (22). Then, Excel spreadsheets can be used for calculations and drawings. For example, Figure 10, Figure 11 and Figure 12 below show Combination (7) in Case 2 of the previous section.

Figure 10, Figure 11 and Figure 12 show that according to the explicit solution of Equation (34), the simulation results obtained by using Excel are consistent with the results obtained by applying MATLAB in the previous section. Applying Excel calculations and graphing may be easier for the water authority to accept than MATLAB.

7. Conclusions

The control design of urban comprehensive water consumption systems with the combination of an urban water demand function based on per capita water consumption and the C-D production function can serve as a reference for optimizing water conservancy input.

The comparison between the model in this paper and some economic dynamic models, such as the Solow model, the Shell model and the aggregate model with renewable labor resources, which are used in the economic, ecological and environmental sciences, showed that the state equations of capital stock and labor stock were different from the above existing models. This eliminated value constraints of these models on investment and labor input, and avoided situations where these models are difficult to apply to urban water supply agencies.

In contrast to the optimal control method, the feedback linearization control method can obtain an explicit solution of the control law, but does not need to solve the complex two-point boundary value problems of differential equations, making this method easier than the optimal control method in application.

Two different cases of urban comprehensive water consumption systems, in the growth period and the declining period, were simulated. The simulation results verified the validity of the proposed control scheme.

In the case simulation of this paper, the iterative solution using the fourth-order and fifth-order Runge-Kutta-Fehlberg algorithm with a variable step size was calculated in MATLAB. More intuitively, it can also be calculated in an Excel spreadsheet in accordance with the corresponding explicit solution. This method is more direct and easier to operate, and perhaps more desirable to a governmental planning committee.

Finally, it should be pointed out that there are also some limitations in this study. First, the process and speed of system stability is dependent on the selection of values for $k_1$ and $k_2$. When the values of other parameters like $D*$ and $S*$ of the system in Equation (22) vary, there is no simpler method to directly determine the appropriate values of $k_1$ and $k_2$. Instead, we need to combine the different values of $k_1$ and $k_2$, verify the effects of these combinations on the system operation through simulation, and then determine the values of $k_1$ and $k_2$, which consumes more time.

Author Contributions

K.L. conceived and designed the model; K.L. and T.M. performed the case simulation and analyzed the data; K.L. wrote the paper; T.D. derived and calculated the explicit solution, and T.D. and G. W. reviewed and edited the manuscript. All coauthors contributed substantially to the work reported.

Funding

This work was supported by the National Social Science Fund of China (grant number 17BGL220).

Conflicts of Interest

The authors declare no conflicts of interest.

Figures

Figure 1. Simulation results of Combination (1).

Figure 2. Simulation results of Combination (2).

Figure 3. Simulation results of Combination (3).

Figure 4. Simulation results of Combination (4).

Figure 5. Simulation results of Combination (5).

Figure 6. Simulation results of Combination (6).

Figure 7. Simulation results of Combination (7).

Figure 8. Simulation results of Combination (8).

Figure 9. Simulation results of Combination (9).

Figure 10. D(t) and S(t) in Combination (7) of Case 2.

Figure 11. I(t) in Combination (7) of Case 2.

Figure 12. P(t) in Combination (7) of Case 2.