Model predictive control of agro-hydrological systems based on a two-layer neural network modeling framework

Based on first principles, the soil water dynamics can be modeled using the 1-dimensional Richards equation with the crop and weather information taken into account [11]:

where $h$ $[m]$ is the capillary potential, $c$ $[m^{-1}]$ denotes the soil capillary capacity, and $K$ $[m/s]$ is the soil hydraulic conductivity. The $-\alpha(h)\frac{K_{c}ET_{0}}{z_{r}}$ term captures the effect of weather and crop, where $\alpha(h)$ and $K_{c}$ are dimensionless factors, namely the water stress factor and the crop coefficient. $ET_{0}$ $[m]$ denotes the reference evapotranspiration and is assumed to be a known weather-dependent disturbance. $z_{r}=-0.13$ m represents the rooting depth of the crop and is assumed to be time-invariant in this work. Note that $z$ $[m]$ represents the vertical axis of the soil with upwards as the positive direction. $c(h)$ and $K(h)$ are modeled by the van Genuchten-Mualem soil hydraulic model [19]:

where $K_{s}$ $[m/s]$ denotes the saturated hydraulic conductivity, $\theta_{s}$ $[m^{3}/m^{3}]$ and $\theta_{r}$ $[m^{3}/m^{3}]$ represent the saturated and the residual soil moisture, respectively. $\alpha$ $[m^{-1}]$ and $n$ are the van Genuchten-Mualem parameters related to the soil properties. Sandy loam soil is considered in this work and the corresponding parameters are adopted from [20] and are presented in Table 1.

The boundary conditions employed in this work are presented in (4)–(5). The Neumann boundary condition is used at the soil surface as shown in (4), while free drainage is assumed for the bottom boundary as described in (5).

where $I(t)$ $[m/s]$ and $P(t)$ $[m/s]$ denote the rate of irrigation and precipitation at time $t$. Note that interception caused by grass leaves is ignored in this work, only evapotranspiration, precipitation, and irrigation are considered.

The system of equation (1) is discretized over the vertical spatial axis following the approach proposed in [11], which leads to an ODE system:

where $x\in\mathbb{R}^{26}$ and $u\in\mathbb{R}^{1}$ are the state and input vectors respectively, while $f$ defines the nonlinear system dynamics. The total depth of soil considered is $0.5$ m and is evenly discretized into 26 nodes. The water dynamics at the center of each node are used to represent the dynamics of the entire node. In this work, irrigation rate $I(t)$ $[m/s]$ is the system input, and $h$ $[m]$ at the discretized nodes represents the system states. The system output $y\in\mathbb{R}^{1}$ is the volumetric water content $\theta$ $[m^{3}/m^{3}]$ at the rooting depth $z_{r}$, which is an algebraic function of the state vector defined by the soil-water retention equation of van Genuchten (7). Eqn. (6) and Eqn. (7) are used to generate all training datasets.

In this work, LSTM NNs will be used to approximate the original nonlinear system dynamics due to their ability of efficiently handle sequential data [21, 22]. The motivation for designing the proposed framework is that a single LSTM is not able to capture the complex nonlinear dynamics of the agro-hydrological system over the entire operating range of interest. It was noticed, however, within a smaller operating region, a single LSTM provides satisfying modeling performance. Detailed simulation results will be presented in Section 4.

A schematic diagram of the proposed framework is presented in Figure 2. The first layer contains $K$ sub-models, namely Model $1$– Model $K$, which will be referred to as $M_{1},M_{2},\cdots,M_{K}$ in the following text. The second layer consists of a single model (Model A, referred to as $M_{A}$). Each sub-model in the first layer is responsible for a sub-operating region. The union set of the sub-operating regions forms the entire operating region of interest. $M_{1}-M_{K}$ are LSTMs while $M_{A}$ is a fully connected NN. At the current time step $t$, $M_{1}-M_{K}$ take the system input $u$ and output $y$ obtained from $p$ previous time steps as NN inputs:

and returns $K$ one-step-ahead predictions regarding the system output:

Recall that each sub-model is only reliable in a sub-operating region. Inspired by Mao et al. [8], the proposed framework adopts the idea of the LPV models. The model performance can be improved by combining predictions made by multiple sub-models. The one-step-ahead predictions made by $M_{1}-M_{K}$ are fed to $M_{A}$, where the final prediction of the system output at the next time step ($NN_{out_{2}}=y_{pred}(t+1)$) is obtained.

Overall, at any time step $t$, the proposed framework takes $NN_{in}(t)$ as input and returns the one-step-ahead prediction of the system output $y(t+1)$, which can be expressed as a nonlinear function:

As discussed in the previous section, the first layer of the proposed framework consists of $K$ LSTMs. In this work, $K=3$ is used. The operating regions in terms of the system output for each sub-model are summarized in Table 2. The number of LSTM layers employed in the models is also included. Each model consists of 1 or 2 LSTM layers, followed by a fully connected layer. Activation functions ‘sigmoid’ and ‘tanh’ are used. Both activation functions have smooth dynamics, where ‘tanh’ function provides steeper responses compare to ‘sigmoid’ function. Each sub-model is equivalent to a nonlinear function:

where the model input is a matrix with $20$ rows and $2$ columns, and the output is a scalar. $20$ is the number of the past input-output data points ($p=20$), while $2$ is the number of system input ($n=1$) plus the number of system output($l=1$). All models have a sampling time of 2 hours.

Three training datasets that vary in the corresponding operating region are generated for the sub-models respectively. Each dataset contains 30000 data points. Multi-level pseudorandom input signals are employed to generate all three datasets, which are fed into the discrete ODE system (6) and the algebraic equation (7) for open-loop simulations. A multi-level pseudorandom signal is similar to a pseudorandom binary signal, except that it has multiple levels, which helps to stimulate the nonlinear dynamics of the agro-hydrological system. More details regarding the design of multi-level pseudorandom inputs can be found in our previous work [18]. Figure 3 presents a segment of the training input signal and the corresponding system output for $M_{2}$. Note that a random noise signal $\epsilon$ is added to the output trajectory to mimic real process operations, which is shown in Figure 3(b) in red. The max magnitude of $\epsilon$ is defined as follows:

where $y_{max}$ and $y_{min}$ are the maximum and minimum value of the system output in the dataset.

All datasets need to be reorganized to match the input-output structure of the NNs. Scaling is also required such that the magnitudes of the scaled data are less than 1. All simulations are carried out in Python. TensorFlow [23] is used to train the NN models. Readers may refer to our previous work [18] for more details regarding data pre-processing and NN training.

The second layer of the proposed framework consists of a fully connected NN equivalent to the following nonlinear function:

The inputs to model $M_{A}$ are the one-step-ahead predictions of the system output made by $M_{1}$, $M_{2}$ and $M_{3}$, while the output is the actual system output prediction. Two fully connected layers with the activation function ‘sigmoid’ are used in $M_{A}$.

Under the objective of predicting the system output over the entire operating range of interest, the training dataset of $M_{A}$ covers the entire operating range. A section of the training dataset including 100000 points is presented in Figure 4, where the input signal and the corresponding output trajectories without and with white noise added are presented. Note that the training dataset of $M_{A}$ spans all the sub-operating regions with more aggressive dynamics. Instead of continuous irrigation with lower magnitudes, an impulse irrigation signal is used to provide a more aggressive output response that varies in a larger operating range, which is desired for $M_{A}$ training.

Similar to that discussed in the previous section, reconstruction and scaling are necessary for this dataset. Note that the inputs to $M_{A}$ are not the system input and output. The dataset is passed through the pre-identified sub-models after the data pre-processing. The predictions made by the sub-models are collected, which are the designed inputs for $M_{A}$ and are used to train $M_{A}$.

The multi-step-ahead prediction performances of $M_{1}-M_{3}$ are presented in this section. Validation datasets are obtained in similar manners as the training datasets with 20% white noise added to the system output. This helps to better test the prediction performance of the models. At any given time step, noise-treated validation data are fed to the model for one-step-ahead prediction. Starting from the second time instant, the prediction made by the model at the previous time instant will be used as the initial condition of the current time instant. This procedure is repeated 20 times, providing twenty-steps-ahead predictions.

The prediction performance is measured by the normalized root mean square error (NRMSE) of the multi-step-ahead predictions in the desired operating range of each model and is summarized in Table 3. The NRMSE is defined as follows:

where $t_{f}$ is the total number of validation data points, $y_{act}$ and $y_{pred}$ are the true and predicted system outputs, respectively. Figures 5, 6, and 7 present the multi-step-ahead open-loop prediction performance of $M_{1}$, $M_{2}$, and $M_{3}$ respectively. The same line styles are used in Figures 5–7 for consistency.

All three models provide reasonable multi-step-ahead prediction without delays in the presence of noise. When the magnitude of the system output is low, small prediction offsets can be observed for $M_{1}$ in Figure 5. Figure 6 shows only some minor offsets exist at the local extremes for $M_{2}$. For both $M_{1}$ and $M_{2}$, the effect of noise seems to be minor. For $M_{3}$, Figure 7 indicates that the noise has a stronger impact on the single-step-ahead prediction, which deviates and becomes insignificant as the number of prediction steps increases. It is also noticed that when the output converges to steady-states, the impact of the noise is the most significant. Similar to $M_{1}$ and $M_{2}$, $M_{3}$ also has some minor offsets at the local extremes.

For different operating ranges, different system output dynamics are observed. When the soil moisture content and the irrigation rate are low, no obvious steady-state is observed (Figure 5). As the operating range shifts to the higher magnitude side, the soil moisture content tends to saturate to steady-states under continuous irrigation (Figure 7).

The multi-step-ahead prediction performance of the proposed framework is investigated in this section. A single LSTM that aims to predict the agro-hydrological dynamics over the entire operating region is used as the benchmark approximation model. The benchmark LSTM has the same input ($NN_{in}$) and output ($NN_{out_{2}}$) structures as the proposed two-layer NN framework, which are discussed in Section 3.1. The single LSTM has 2 LSTM layers followed by a fully connected output layer. Both the LSTM layers use the activation function ‘sigmoid’ while the output layer uses the activation function ‘tanh’. Note that the dataset employed for training the benchmark LSTM is the same as that used to train the two-layer NN model.

The same validation dataset is used to compare the performance of the two approximation approaches, which spans the entire operating range of interest. The open-loop prediction trajectories of the system output based on the benchmark LSTM and the proposed framework are presented in Figures 8 and 9 respectively. The single-step-ahead, ten-step-ahead, and twenty-step-ahead predictions are presented in each figure with the same line styles as those employed in the previous section. The NRMSE values under the two cases and the percentage difference $\triangle$ between them are summarized in Table 4. The percentage difference is calculated with respect to the NRMSE of the benchmark LSTM:

As presented in Figure 8, the LSTM can respond to aggressive dynamics efficiently without delays. However, significant offsets are observed at the local extremes with larger magnitudes. A section of the trajectories is amplified, which can be seen on the top-right corner of Figure 8. At both the local maximum and the minimum, significant offsets can be observed in the zoomed-in plot, of which the predicted local minimums are at the level of the true local maximum. Better prediction performance is observed in Figure 9. To be specific, the dynamics at the local minimums can be captured by the proposed two-layer framework much better compared to the LSTM. From the zoomed-in window on the top right, it can be observed that the predictions of the local minimum and maximum are more accurate compared to the LSTM. Table 4 indicates comparison results, showing that the NRMSE values of the proposed framework are significantly lower than those obtained with the LSTM for all different prediction steps. The improvement of the proposed framework is the most significant with one-step-ahead predictions and slowly reduces as the prediction steps increase.

It is noticed that the proposed framework cannot accurately approximate the local maximums with larger magnitudes. To account for this problem, we introduce a bias-and-scale-updated model to capture the plant-model-mismatch, of which the details are provided in the following section.

In this section the performance of the proposed framework in ZMPC applications is validated by comparing the simulation results with those obtained based on the discretized Richards equation and the single LSTM. Table 6 summarizes the average computational time for solving the optimization problem, the economic performance of ZMPC simulations, and the percentage difference in economic performance based on different models. The economic performance is measured by the accumulated irrigation amount $I_{T}$ in millimeters, which is calculated as follows:

where $u(t)$ is the irrigation rate at time $t$ with the unit $m/s$. The percentage difference in the accumulated irrigation is calculated with respect to the result obtained based on the Richards equation. The optimal control strategy and the corresponding output trajectories obtained based on different models are presented in Figure 10. The results shown in Table 6 and Figure 10 are obtained in the presence of 2% process noise. Figure 11 presents the optimal system input and output trajectories obtained with 5% process and 5% measurement noise based on the single LSTM and the proposed framework. The colors and markers employed in Figure 11 are the same as those in Figure 10. Table 7 summarizes more ZMPC simulation results based the single LSTM and the proposed framework in the presence of process and measurement noise.

As presented in Table 6, both the single LSTM and the proposed framework help to reduce the computational cost significantly. The computational cost reduction achieved by the single LSTM is more significant compare to the two-layer NN due to the simplicity of its structure. With only 8.66% deviation, the economic performance of controller using the two-layer NN is similar to that using the Richards equation, which is significantly better compared to that obtained with the single LSTM. The trajectories in Figure 10 indicate the same, as the optimal control law and the system response of the two-layer NN and the Richards equations are highly overlapped. Based on the Richards equation, no irrigation is applied for the first a few steps and starts to increase slowly over time. The two-layer NN leads to similar control strategy with the minor differences caused by the plant-model-mismatch. As for the single LSTM, the irrigation policy is less efficient, which is nonzero over the entire horizon. The output trajectories converge to the lower bound of the target zone for simulations based on the Richards equation and the proposed framework. The responses are reasonable as the economic objective is to reduce the amount of irrigation required. On the other hand, for the single-LSTM-based ZMPC, slower convergence and offsets with respect to the lower bound of the target zone is observed in the optimal output trajectory.

From Tables 6 and 7, it is observed that compared to process noise, measurement noise has stronger impacts on the simulation results. As shown in Table 7, except for the case that considered 5% process noise and 5% measurement noise, control with the two-layer NN has significantly better economic performances compared to that using the single LSTM. As shown in Figure 11, in the presence of 5% process noise and 5% measurement noise, the two-layer NN provides more aggressive control policy to capture the effect of the noise compare to the single LSTM. However, the impact of noise in this case is more significant compare to the control action due to the slow-responding nature of the agro-hydrological system, as the output trajectories obtained based on the two models are similar. To summarize, the proposed framework is shown to have better performances compared to the single LSTM in ZMPC applications. However, some minor violations of the target zone is observed at the end of the horizon of interest (Figure 11) due to the effect of noise, which is not ideal as low soil moisture may cause severe consequences in crop growing. Furthermore, precipitation is an essential disturbance in real-world applications and should be considered. In the following section, we investigate the performance of a shrinking target zone in capturing the effect of significant noise and weather disturbances.

In this section, the performance of the proposed shrinking-zone ZMPC and the ZMPC with time-invariant zone (referred to as the basic ZMPC in the following text) is compared in the presence of 5% measurement noise and weather disturbance. The controller is assumed to have information regarding the weather through weather forecasting with some uncertainties. In this work we consider 20% error in weather forecasting. Various parameter combinations are investigated. In Section 7.2.1, $\mu$ is kept constant and the effect of $\underline{Y}$ and $\overline{Y}$ discussed. The opposite is done in Section 7.2.2, where $\underline{Y}$ and $\overline{Y}$ are kept constant with varying $\mu$.