Grey-Box Modeling and Decoupling Control of a Lab Setup of the Quadruple-Tank System

: The quadruple-tank system (QTS) is a popular educational resource in universities for studying multivariable control systems. It enables the analysis of the interaction between variables and the limitations imposed by multivariable non-minimum phase zeros, as well as the evaluation of new multivariable control methodologies. The works utilizing this system present a theoretical model that may be too idealistic and based on erroneous assumptions in real-world implementations, such as the linear behavior of the actuators. In other cases, an identified linear model is directly provided. This study outlines the practical grey-box modeling procedure conducted for the QTS at the University of Cordoba and provides guidance for its implementation. A configurable nonlinear model was developed and controlled in a closed loop using different controllers. Specifically, decentralized control, static decoupling control, and simplified decoupling control were compared. The simulation designs were experimentally validated with high accuracy, demonstrating that the conclusions reached with the developed model can be extrapolated to the real system. The comparison of these three control designs illustrates the advantages and disadvantages of decoupling in certain situations, especially in the presence of non-minimum phase zeros.


Introduction
Industrial processes are often multivariable and exhibit nonlinear behavior [1].The interactions between components, variables to be controlled, and manipulated variables can increase the complexity of a nonlinear system.Due to the characteristics of nonlinearity, time variation, and cross-coupling between variables, classical control methods developed for single-variable systems may lead to inefficient results when they are applied to multivariable processes [2].Additionally, the design of a multivariable process control system must consider the possible presence, location, and direction of multivariable zeros, as they can affect system stability.The design of control systems is generally based on a process model in an iterative procedure to save time and money.It is important to maintain a clear and logical flow of information with causal connections between statements.Although closed-loop control systems are less sensitive to process modeling errors, obtaining a proper model is crucial for designing a good controller.This is especially true for multivariable processes, which present greater control difficulties that the model should reflect.
Industrial process modeling has traditionally been based on white-box modeling or black-box identification.White-box modeling involves building a model using scientific relationships that fully describe the process.On the other hand, black-box identification employs a parametric model that fits the measured process data obtained experimentally.For many industrial processes, there is some knowledge available, although incomplete, about the system.Between white-box and black-box models, there exists a grey area that offers a third way to create engineering system models.This grey-box approach utilizes a priori knowledge about the process and estimates the unknown parts of the model from experimentally measured data [3].When there is good prior knowledge of the mechanisms underlying the behavior of a process, the relevant equilibrium equations can be expressed as a set of differential equations.In the simplest case, the system may have some unknown parameters that need to be estimated from measured data.In the general case, the model's initial structure is formulated on the basis of physical knowledge and then refined to fit the experimental data [4].The process of developing a mechanistic grey-box model is typically oriented toward objects or sub-models and involves a series of steps that require solving less-demanding subproblems, such as basic modeling, experiments, estimation, data fitting, and model validation.Thus, grey-box models are simplified representations of the different components that facilitate system integration [5].
After obtaining the process model, and given the control specifications, it is necessary to decide which control structure to use.For multivariable processes, there are typically two approaches: decentralized control and centralized control [6].The interaction of a multiple-input multiple-output (MIMO) process occurs when one input variable affects several output variables to a greater or lesser degree, making it challenging to operate the process and design its control system.Therefore, any multivariable control design methodology should consider interaction.When the interaction is insignificant, control can be approached using a decentralized (also known as multi-loop) scheme.This scheme uses independent controllers as if the process were composed of monovariable processes [7], as shown in Figure 1a.The references signals are represented by r i , the controlled variables by y i , and the control signals by u i .These controllers may require slight retuning of their gains.In other cases where there is more interaction, and this is taken into account, the decentralized structure is still used.However, each controller is designed to incorporate the effects of the other loops on the corresponding process [8].Before designing a decentralized control, it is important to choose a suitable pairing between input and output variables that presents as little interaction as possible.The relative gain array (RGA) [9] is a commonly used measure to solve pairing problems.When the interaction is significant, centralized control is necessary, as shown in Figure 1b.This involves designing a single control block that takes into account all the information from the measured variables to generate all the control signals.However, a common intermediate solution in the industry is to incorporate a decoupling network (or decoupler) in series with the process [6].The intention is to eliminate or reduce interaction and indirectly facilitate the design of a decentralized control system (see Figure 1c).
For many industrial processes, there is some knowledge available, although incomplete, about the system.Between white-box and black-box models, there exists a grey area that offers a third way to create engineering system models.This grey-box approach utilizes a priori knowledge about the process and estimates the unknown parts of the model from experimentally measured data [3].When there is good prior knowledge of the mechanisms underlying the behavior of a process, the relevant equilibrium equations can be expressed as a set of differential equations.In the simplest case, the system may have some unknown parameters that need to be estimated from measured data.In the general case, the model's initial structure is formulated on the basis of physical knowledge and then refined to fit the experimental data [4].The process of developing a mechanistic grey-box model is typically oriented toward objects or sub-models and involves a series of steps that require solving less-demanding subproblems, such as basic modeling, experiments, estimation, data fitting, and model validation.Thus, grey-box models are simplified representations of the different components that facilitate system integration [5].
After obtaining the process model, and given the control specifications, it is necessary to decide which control structure to use.For multivariable processes, there are typically two approaches: decentralized control and centralized control [6].The interaction of a multiple-input multiple-output (MIMO) process occurs when one input variable affects several output variables to a greater or lesser degree, making it challenging to operate the process and design its control system.Therefore, any multivariable control design methodology should consider interaction.When the interaction is insignificant, control can be approached using a decentralized (also known as multi-loop) scheme.This scheme uses independent controllers as if the process were composed of monovariable processes [7], as shown in Figure 1a.The references signals are represented by ri, the controlled variables by yi, and the control signals by ui.These controllers may require slight retuning of their gains.In other cases where there is more interaction, and this is taken into account, the decentralized structure is still used.However, each controller is designed to incorporate the effects of the other loops on the corresponding process [8].Before designing a decentralized control, it is important to choose a suitable pairing between input and output variables that presents as little interaction as possible.The relative gain array (RGA) [9] is a commonly used measure to solve pairing problems.When the interaction is significant, centralized control is necessary, as shown in Figure 1b.This involves designing a single control block that takes into account all the information from the measured variables to generate all the control signals.However, a common intermediate solution in the industry is to incorporate a decoupling network (or decoupler) in series with the process [6].The intention is to eliminate or reduce interaction and indirectly facilitate the design of a decentralized control system (see Figure 1c).The decoupling of variables, or the absence of interaction, is a desirable feature in many industrial applications.This makes it easier for technicians to determine the reference values needed to achieve specific objectives in a multivariable control system.Additionally, it allows for independent improvement of system responses.Leading manufacturers of distributed control systems [10] consider interaction one of the main problems in industrial processes.Historically, they have supported decoupling techniques as a means of combating it.In [11], decoupling is considered an advanced regulatory control element that can be used to improve control in cases where a simple control loop is not sufficient.
Traditionally, control formulations for decoupling have focused on 2 × 2 processes, which have two inputs and two outputs.This is because these processes are either the most common or because more complex processes are often decomposed into 2 × 2 blocks with important interactions between their inputs and outputs.Decoupling control assumes that the process model is represented by a nonsingular square matrix of transfer functions G(s), with no poles in the right half-plane (RHP).The decoupler is represented by a square transfer matrix D(s), and the decentralized controller is represented by a diagonal matrix C(s). Figure 2 illustrates the general scheme of a 2 × 2 control system that combines a decoupler and a diagonal controller.The decoupling network is placed between the diagonal controller and the process to form the new apparent process, Q(s) = G(s)D(s).The diagonal controller, C(s), manipulates the variables v i instead of the real manipulated variables u i .The compensating block D(s), also known as the decoupler, is designed to eliminate or reduce the interaction between the variables of the apparent process.This allows the controller C(s) to perceive it as a set of totally independent processes or with significantly less interaction.Monovariable controllers c i (s) would be designed for each of these new apparent processes to form the diagonal control C(s) [6,10,12].The decoupling of variables, or the absence of interaction, is a desirable feature in many industrial applications.This makes it easier for technicians to determine the reference values needed to achieve specific objectives in a multivariable control system.Additionally, it allows for independent improvement of system responses.Leading manufacturers of distributed control systems [10] consider interaction one of the main problems in industrial processes.Historically, they have supported decoupling techniques as a means of combating it.In [11], decoupling is considered an advanced regulatory control element that can be used to improve control in cases where a simple control loop is not sufficient.
Traditionally, control formulations for decoupling have focused on 2 × 2 processes, which have two inputs and two outputs.This is because these processes are either the most common or because more complex processes are often decomposed into 2 × 2 blocks with important interactions between their inputs and outputs.Decoupling control assumes that the process model is represented by a nonsingular square matrix of transfer functions G(s), with no poles in the right half-plane (RHP).The decoupler is represented by a square transfer matrix D(s), and the decentralized controller is represented by a diagonal matrix C(s). Figure 2 illustrates the general scheme of a 2 × 2 control system that combines a decoupler and a diagonal controller.The decoupling network is placed between the diagonal controller and the process to form the new apparent process, Q(s) = G(s)D(s).The diagonal controller, C(s), manipulates the variables vi instead of the real manipulated variables ui.The compensating block D(s), also known as the decoupler, is designed to eliminate or reduce the interaction between the variables of the apparent process.This allows the controller C(s) to perceive it as a set of totally independent processes or with significantly less interaction.Monovariable controllers ci(s) would be designed for each of these new apparent processes to form the diagonal control C(s) [6,10,12].The decoupler design can be approached as a straightforward mathematical problem.For the 2 × 2 case, perfect direct decoupling can be achieved by formulating it as shown in (1), where the apparent process Q(s) should only have elements on the main diagonal; the Laplace variable s has been omitted.
When designing a decoupling network, there are certain degrees of freedom for choosing the decoupling elements, resulting in different decouplers, some simpler and others more complex to implement [13].One of the most common decouplers is the simplified decoupling [14], which involves fixing two elements of D(s) to unity, usually the diagonal ones.The resulting decoupler is shown in (2), along with the expressions of the resulting apparent processes Q(s).In any case, the decoupling elements must be realizable; therefore, they must be proper, causal, and stable.Decoupling control design approaches are based on a process model and may have robustness issues in certain nonlinear The decoupler design can be approached as a straightforward mathematical problem.For the 2 × 2 case, perfect direct decoupling can be achieved by formulating it as shown in (1), where the apparent process Q(s) should only have elements on the main diagonal; the Laplace variable s has been omitted.
When designing a decoupling network, there are certain degrees of freedom for choosing the decoupling elements, resulting in different decouplers, some simpler and others more complex to implement [13].One of the most common decouplers is the simplified decoupling [14], which involves fixing two elements of D(s) to unity, usually the diagonal ones.The resulting decoupler is shown in (2), along with the expressions of the resulting apparent processes Q(s).In any case, the decoupling elements must be realizable; therefore, they must be proper, causal, and stable.Decoupling control design approaches are based on a process model and may have robustness issues in certain nonlinear systems.This may require a redesign of the control system if the operating point changes.Recent works propose model-free robust U-control approaches for decoupling [15] and decentralized control [16]  (2) The quadruple-tank system (QTS) is a popular choice for teaching multivariable system concepts and testing different multivariable control methodologies, including decoupling.It was proposed by [17][18][19] and is widely used in university laboratories.The system involves controlling the water level in tanks that are coupled together.Controlling liquid levels in systems with interconnected tanks can be a challenging multivariable control problem.This issue is relevant in various industries, including water treatment, chemical, and biochemical plants, food processing, metallurgy, filtration devices, etc. [2,20].The QTS is a nonlinear system that illustrates several phenomena of multivariable systems, such as the interaction and effects of multivariable zeros [21].Therefore, it presents itself as an excellent tool for teaching multivariable control techniques and understanding their advantages and disadvantages.The QTS has been used to implement multiple control methods, both classical and advanced, for teaching and research purposes.In one of the first studies on QTS [19], a decentralized PI control was designed for the QTS configured with a multivariable RHP zero.Other authors have applied internal model control [22], multivariable H ∝ control [3], quantitative feedback control [23], LQG optimal control [24], predictive control [25,26], and distributed model predictive control [27].More recent works have applied nonlinear techniques to the QTS such as sliding mode control [28,29], feedback linearization [20], fuzzy control [30,31], and neural networks [32], among others.
In all of these studies, the presented theoretical model of the QTS is simple and intuitive.However, experience has shown that it is based on certain assumptions that may be incorrect in real implementations of this process.In most cases, the nonlinear model presented is very idealistic, assuming a linear behavior of the system actuators, such as three-way valves or pumps.In addition, little information is provided on how the model parameters are obtained.In some cases, the QTS is presented directly as a linear model obtained through identification at an operating point.This work presents a detailed description of the practical grey-box modeling procedure for the QTS at the University of Cordoba, including practical aspects that have not been previously described in other works.This procedure can also be applied to other systems with liquid tanks.Additionally, suggestions are provided for constructing the QTS.A configurable nonlinear model with good accuracy is obtained, allowing for the adjustment of the system to show little or large interactions or the position of a multivariable zero.The model is validated in a closed loop with different controllers.Specifically, decentralized PI control, static decoupling control, and simplified decoupling control are compared.The simulation designs were experimentally validated with high accuracy.Furthermore, the comparison of the three control designs demonstrates the benefits and drawbacks of decoupling under specific circumstances.
The rest of the paper is structured as follows: In Section 2, the experimental quadrupletank system located in the Process Control Laboratory of the University of Cordoba is described.The system was implemented to evaluate different multivariable control methodologies in a real process.This section covers the operation of the system, as well as the components and instrumentation used in its construction.Section 3 determines the nonlinear model of the plant through grey-box modeling.The resulting expressions of the linear model, linearized around a generic operating point, are derived.Section 4 illustrates key aspects of multivariable control through the presentation of the three proposed control strategies applied to QTS.Section 5 discusses the simulation and experimental results of the three control systems.Finally, Section 6 summarizes the main conclusions.

Laboratory Setup of the Quadruple-Tank System
The quadruple-tank system implemented in the Process Control laboratories of the University of Cordoba follows the scheme depicted in Figure 3a.The system consists of a reservoir tank at the bottom, and four tanks located on two levels: tanks 1 and 2 on the lower level, and tanks 3 and 4 on the upper level.The pump on the left (p L ) delivers a flow of water f L from the reservoir tank, which is then split into two flows through the three-way valve (V L ): f 1 to tank 1 and f 4 to tank 4. The pump on the right (p R ) drives a flow rate f R , which is then split by the V R valve into f 2 , which goes to tank 2, and f 3 , which goes to tank 3. The tanks drain through manual valves that can be multiturn (v 1 , v 2 , v 3 , and v 4 ) or on/off (v 1B , v 2B , v 3B , and v 4B ).Tanks 1 and 2 divert water to the reservoir tank, whereas tanks 3 and 4 divert water to tanks 1 and 2, respectively.The tanks are coupled, meaning that, for example, when the pump on the right is used to fill tank 2, a portion of the total flow f R goes to tank 3 and eventually reaches tank 1.The variables measured by the sensors are the water levels in the tanks (h 1 , h 2 , h 3 , and h 4 ) and the total flow rates of each branch (f L and f R ).The plant actuators are the frequency inverters that drive the pumps; therefore, the variables that can be manipulated are the frequency setpoints of these inverters, expressed as a percentage from 0 to 100%.However, as explained below, a secondary flow control loop is closed in each branch (left and right), so that the new manipulated signals are the flow references of each branch f L_ref and f R_ref .
Actuators 2024, 13, x FOR PEER REVIEW 5 of 22 control strategies applied to QTS.Section 5 discusses the simulation and experimental results of the three control systems.Finally, Section 6 summarizes the main conclusions.

Laboratory Setup of the Quadruple-Tank System
The quadruple-tank system implemented in the Process Control laboratories of the University of Cordoba follows the scheme depicted in Figure 3a.The system consists of a reservoir tank at the bottom, and four tanks located on two levels: tanks 1 and 2 on the lower level, and tanks 3 and 4 on the upper level.The pump on the left (pL) delivers a flow of water fL from the reservoir tank, which is then split into two flows through the threeway valve (VL): f1 to tank 1 and f4 to tank 4. The pump on the right (pR) drives a flow rate fR, which is then split by the VR valve into f2, which goes to tank 2, and f3, which goes to tank 3. The tanks drain through manual valves that can be multiturn (v1, v2, v3, and v4) or on/off (v1B, v2B, v3B, and v4B).Tanks 1 and 2 divert water to the reservoir tank, whereas tanks 3 and 4 divert water to tanks 1 and 2, respectively.The tanks are coupled, meaning that, for example, when the pump on the right is used to fill tank 2, a portion of the total flow fR goes to tank 3 and eventually reaches tank 1.The variables measured by the sensors are the water levels in the tanks (h1, h2, h3, and h4) and the total flow rates of each branch (fL and fR).The plant actuators are the frequency inverters that drive the pumps; therefore, the variables that can be manipulated are the frequency setpoints of these inverters, expressed as a percentage from 0 to 100%.However, as explained below, a secondary flow control loop is closed in each branch (left and right), so that the new manipulated signals are the flow references of each branch fL_ref and fR_ref.Figure 3b shows a photograph of the real system.Its construction, unlike other similar plants, prioritized the use of industrial components over esthetics, resulting in larger dimensions.The elements that comprise it can be grouped into two categories: the Figure 3b shows a photograph of the real system.Its construction, unlike other similar plants, prioritized the use of industrial components over esthetics, resulting in larger dimensions.The elements that comprise it can be grouped into two categories: the hydraulic circuit and control system [33].The hydraulic circuit is supported by a metallic structure and is composed of the following: • Two 0.30 kW three-phase centrifugal pumps.

•
Five methacrylate tanks.The reservoir tank measures 60 × 50 × 40 cm.The other four tanks have a square internal section of 19 × 19 cm 2 and a height of 35 cm for the lower tanks and 70 cm for the upper tanks.To improve the response time of the tanks, cylindrical tubes of constant section were introduced, as shown in Figure 3b.The cross-sectional area of the lower tanks, where only one tube is present, is 316.82 cm 2 .The cross-sectional area of the upper tanks, where four tubes are present, is 184.29 cm 2 .• Two manual three-way diverter valves (V L and V R ) are used to adjust the proportion of flow going to the lower tank in each branch.When in the 0% position, all the flow goes to the corresponding upper tank, so no water reaches the lower tank of that branch.Conversely, when in the 100% position, all the flow goes to the lower tank.

•
Eight manual two-way valves are used to adjust the outflows of the four tanks.The multi-turn valves v 1 , v 2 , v 3 , and v 4 can vary their opening degree with a total of up to six and a half or seven turns, depending on the valve.In contrast, valves v 1B , v 2B , v 3B, and v 4B have only two positions: fully closed or fully open.• We use 10/8 mm polyurethane pipes and pipe fittings for connecting the different elements of the hydraulic circuit.In addition, two 40/32 mm PVC hoses are used to drain the inlet flows from the three-way valves to the lower tanks.Unlike the schematic shown in Figure 3a, the two outlets from each three-way valve are brought to the same height through 10/8 mm pipes.One branch goes to the corresponding upper tank, whereas the other branch goes to the 44/38 mm pipe so that from the upper height it drains without head loss.This design ensures that the two outlet branches of the three-way valves have similar head loss and that the flow distribution, depending on the position of the valve, is independent of the total inlet flow to the three-way valve.
The control system comprises measurement sensors, actuation devices, and control elements.These are as follows:

•
Four upwelling diaphragm pressure transmitters for general applications.A measuring range of [0-0.1]bar, which is approximately a [0-100] cm water column.These transmitters were used to measure the water level in the four tanks.• Two electromagnetic flow meters for measuring the total flow rates of the branches f L and f R .The measurement range is set to [0-200] cm 3 /s.• Two frequency inverters to control the flow rate delivered by the pumps through a secondary control loop.• An NI-DAQ 6035E data acquisition card to connect the process signals to a computer.
• A personal computer contains the acquisition card and the control software, which in this case is MATLAB 7.4 and its Real Time Windows Target toolbox.The configured sampling time is 1 s.

Modeling of the Quadruple-Tank System
This section explains the practical procedure proposed for obtaining a nonlinear model of the experimental plant from the modeling of its components using a grey-box approach.Two parts are considered: 1.
Pump dynamics and flow control: they determine the total flow rates of each branch (left and right) that the pumps drive from the reservoir tank to the other tanks.

2.
Dynamics of the four tanks: They determine the water level of each tank as a function of its corresponding inlet and outlet flow rates.The inlet flow rates depend on the total flow rate of each branch, the position of the three-way valves, and in the case of the lower tanks, the outlet flow rates of the upper tanks.The outflow of each tank depends on its water level and the degree of opening of its outlet valves.
Once the complete nonlinear model of the system is obtained and a configuration with the valves is set, it is linearized around an operating point to obtain a linear model.This linear model is represented as a matrix of transfer functions relating the water levels in tanks 1 and 2 to the total flows delivered by the pumps.

Dynamics of the Pumps and Secondary Flow Control Loops
The flow rates f L and f R sent to each branch (left and right) by the pumps are determined using this dynamic as a function of the percentage of operation of the frequency inverter of the corresponding pump, p L or p R .The transfer function in (3) models the relationship between the flow rate (cm 3 /s) of a branch, F(s), and the percentage of operation (%) of the corresponding pump drive, P(s), at approximately 50% operation.
This transfer function may not be entirely accurate because of its nonlinear flow-pump characteristic, which varies depending on the operating point.However, its dynamics are fast enough compared to those of the tanks.Therefore, a secondary flow control loop was implemented in each branch using a proportional-integral (PI) controller to track flow references and reject disturbances.In addition, this approach prevents issues with pump performance fluctuation caused by changes in the reservoir tank level or pump heating.The implemented PI controller has a proportional gain K P of 0.24%/(cm 3 /s) and an integral time constant T I of 2.4 s.These values enable achieving responses with zero position error and a settling time of approximately 15 s, which is sufficiently fast to disregard these dynamics in comparison to those of the tanks.The controller was implemented in MATLAB/Simulink.The new process control signals are the flow references of each branch (f L_ref and f R_ref ) in the range of [0-200] cm 3 /s.These signals are provided from a higher control level.

Tank Dynamics
Tank dynamics determine the inlet and outlet flow rates and the level of each tank.The dynamics are based on Bernoulli's law and the laws of mass balance.The basic dynamics of tank k are given using the differential Equation (4).To determine the evolution of the water level h k , it is only necessary to calculate the inflow f k_in and outflow f k_out , given the area A k of the tank, which in this case is constant. .
The inlet flow rates entering each tank are dependent on the total flow rate and the position of the three-way valve of the corresponding branch.The proportion of flow received by each of the two tanks associated with that branch was experimentally determined for different flow rates and positions of the three-way valve.For this purpose, several experiments were conducted, which consisted of filling the tanks with fixed values of both the reference flow rate and the position of the three-way valve.After the initial transient, the slope of the filling curve was determined to be linear.Multiplying the slope by the known area of the tanks gave the inlet flow rate for each tank.Because the total flow per branch is also known, the fraction going to the lower tanks can be determined.Table 1 displays the ratio γ L (in %) of the flow rate f L directed to tank 1 for varying values of f L and different positions of the left three-way valve V L .Similarly, Table 2 presents the ratio γ R of the flow rate f R directed to tank 2 for different values of f R and various positions of the right three-way valve V R .Figure 4 displays the percentages of these tables in graph form, exhibiting an almost linear trend.
From these tables, the surfaces in Figure 4 are fitted by planes whose equations are given in (5).These equations can be used to determine the fraction of flow that will enter the lower tank for any opening of the three-way valves (ranging from 0 to 100%) and any value of total flow of a branch (in the range 0 to 200 cm 3 /s); the remaining flow will enter the upper cross tank.Thus, using Equation ( 6) and the flow rates of each branch, we can calculate the inlet flow rate for each tank.
(5) value of total flow of a branch (in the range 0 to 200 cm 3 /s); the remaining flow will enter the upper cross tank.Thus, using Equation ( 6) and the flow rates of each branch, we can calculate the inlet flow rate for each tank.
(1 ) The outflow rate of tank k is dependent on the degree of opening of its multiturn valve vk, whether its valve vkB is open or closed, and the water level in the tank.The outflow rate fvk of each valve vk was modeled according to Equation (7), where parameters αk and The outflow rate of tank k is dependent on the degree of opening of its multiturn valve v k , whether its valve v kB is open or closed, and the water level in the tank.The outflow rate f vk of each valve v k was modeled according to Equation (7), where parameters α k and β k vary as a function of the manual opening of the valves, i.e., the number of opening turns.The water level in the tanks is measured from their base.However, because of the position of the tank outlet ports and the pressure sensor, not all of the measurement range is useful.
If the level falls below the minimum value of 2 cm (h min ), the outflow becomes zero or negligible.Additionally, the equation becomes invalid if the tank overflows.
To obtain the coefficients α k and β k for each valve, a series of experiments were performed.The experiments consisted of emptying each tank from its maximum height h k_max for different opening turns of the valve to be modeled.Once the level data in the tank were recorded, the emptying curve h(t) was fitted to a second-degree equation using least squares, as shown in (8).As shown in (9), this equation is derived with respect to time, and the time variable is then solved and substituted into (8).Next, the derivative of h is solved, resulting in an expression that relates the derivative of the tank level to the water level in the tank.Figure 5 shows the real emptying curves (in solid line) and their adjustments (in dashed line) for valves v 1 , v 2 , v 3 , and v 4 when three turns have been opened. .
Actuators 2024, 13, x FOR PEER REVIEW 9 of 22 βk vary as a function of the manual opening of the valves, i.e., the number of opening turns.The water level in the tanks is measured from their base.However, because of the position of the tank outlet ports and the pressure sensor, not all of the measurement range is useful.
If the level falls below the minimum value of 2 cm (hmin), the outflow becomes zero or negligible.Additionally, the equation becomes invalid if the tank overflows.
( ) To obtain the coefficients αk and βk for each valve, a series of experiments were performed.The experiments consisted of emptying each tank from its maximum height hk_max for different opening turns of the valve to be modeled.Once the level data in the tank were recorded, the emptying curve h(t) was fitted to a second-degree equation using least squares, as shown in (8).As shown in (9), this equation is derived with respect to time, and the time variable is then solved and substituted into (8).Next, the derivative of h is solved, resulting in an expression that relates the derivative of the tank level to the water level in the tank.Figure 5 shows the real emptying curves (in solid line) and their adjustments (in dashed line) for valves v1, v2, v3, and v4 when three turns have been opened.As there is no inflow, the outflow of the valve is equal in the absolute value to the derivative of the water volume in the tank.Since the cross-sectional area (A) of each tank is constant, the derivative of the volume and the derivative of the water level are As there is no inflow, the outflow of the valve is equal in the absolute value to the derivative of the water volume in the tank.Since the cross-sectional area (A) of each tank is constant, the derivative of the volume and the derivative of the water level are proportional.Therefore, given the derivative of the tank level, the outflow rate can be determined based on the water height.Thus, the outflow f V can be determined using the expression in (10), which is similar to Equation (7).The parameters α and β are identified as a function of the parameters a, b, and c of the least squares optimization and the tank section A.
Figure 6 displays the estimated values of α and β for the different multiturn valves in the form of circular dots.These values were obtained by conducting the same experiment for various valve openings.For each parameter α k and β k , a polynomial equation was fitted to determine the value as a function of the degree of valve opening.The fitting curves resulting from this analysis are also presented in Figure 6.The polynomial equations are provided in (11), and the fixed coefficients of the on/off valves v 1B , v 2B , v 3B , and v 4B when they are open are given in (12).There is significant symmetry between the valves on the left branch and their corresponding valves on the right side of the system.Most fitting polynomials in (11) are of order 3, but some have a different order to achieve a better fit.For certain parameters, the fitting polynomial of order 3 showed a downward trend at the end points, which does not correspond to reality.As the valve is opened further, both α and β should grow or at least not decrease.proportional.Therefore, given the derivative of the tank level, the outflow rate can be determined based on the water height.Thus, the outflow fV can be determined using the expression in (10), which is similar to Equation ( 7).The parameters α and β are identified as a function of the parameters a, b, and c of the least squares optimization and the tank section A.
Figure 6 displays the estimated values of α and β for the different multiturn valves in the form of circular dots.These values were obtained by conducting the same experiment for various valve openings.For each parameter αk and βk, a polynomial equation was fitted to determine the value as a function of the degree of valve opening.The fitting curves resulting from this analysis are also presented in Figure 6.The polynomial equations are provided in (11), and the fixed coefficients of the on/off valves v1B, v2B, v3B, and v4B when they are open are given in (12).There is significant symmetry between the valves on the left branch and their corresponding valves on the right side of the system.Most fitting polynomials in ( 11) are of order 3, but some have a different order to achieve a better fit.For certain parameters, the fitting polynomial of order 3 showed a downward trend at the end points, which does not correspond to reality.As the valve is opened further, both α and β should grow or at least not decrease.

Complete Nonlinear Model
The water level dynamics in each of the four tanks are calculated based on the differential equation of the tank given in (4), using the inlet and outlet flow rates.The upper tanks receive input only from the pumps, whereas the lower tanks receive input from both the pumps and outflow from their corresponding upper tanks.According to the plant schematic, the set of differential equations is as follows: .

Linearization of the Nonlinear Model
This section derives the linear model in transfer functions by linearizing the nonlinear model around an operating point, based on the modular decomposition of the plant.The general equation in relative variables in the Laplace domain is obtained by linearizing (4) for each tank k around an operating point.
The equation above requires determination of the inlet and outlet flow rates for each tank.These flow rates are the sum of different components, as expressed in (13), depending on the tank.The outlet flow rates for each valve can be linearized using Taylor series expansion and truncated after the first-order term, as expressed in (15).The stationary values of the working point around which the model is obtained are f k_out,0 and h k,0 .After transitioning to the Laplace domain, we obtain Equation ( 16), where f k_out (s) and h k (s) are incremental variables with respect to the operating point.The transfer function is represented by a constant K k_out .
As previously stated, the pump flow rates are determined by secondary control loops that use their flow references as new inputs to the process.The response of these loops is very fast, and they have zero position error due to their integral action; therefore, their dynamics can be ignored.The ratio of flow from each branch to the corresponding lower tank, i.e., the fraction γ L and γ R , is of interest.These ratios depend on the position of the corresponding three-way valve and the total flow rate of its branch according to (5).For example, the partial derivative of γ L with respect to the total flow rate f L_ref of the left branch can be obtained using the expression in (5) for a fixed position of the three-way valve V L .The linearization process is expressed in (17).The transfer function is given by K 1_in , which is another constant.
In the complementary part of that branch, specifically tank 4, the gain would be complementary to one, as shown in (18).Similar models can be obtained in the righthand branch.
To conclude, by linearizing the system in (13), using transfer functions similar to those in ( 14), ( 16)- (18), and solving for h 1 (s) and h 2 (s), we obtain the final expressions of the linear system, given using The transfer matrix of the linearized model can be expressed more concisely as follows: where the outputs are h 1 and h 2 (in cm), and the inputs are f L_ref and f R_ref (in cm 3 /s).The time constants are given using where the general expression of K k_out is given in (16).The steady-state gain matrix and the relative gain matrix (RGA) of the system can be determined using Equations ( 22) and ( 23), respectively.
The zeros of the multivariable process in ( 20) can be determined by finding the zeros of its determinant.According to [16], at least one of these zeros is always located in the left half-plane of the s-plane.If the condition in ( 24) is met, the other zero will be located in the right half-plane (RHP), which can lead to additional control challenges.

Validation of the Model
The experimental plant was configured to have a multivariable RHP zero by adjusting the different valves.Valves v 1 and v 2 are set to be six-turns open, while valves v 3 and v 4 are set to be one-turn open.All on-off valves are open.The three-way valves are set to be in the 30% position to provide more flow to the upper tanks.The reference flow rates selected as inputs were 135 cm 3 /s in both branches.The approximate stationary point is h 1 = 18 cm, and h 2 = 19 cm.The 2 × 2 matrix below represents the linear model obtained from the linearization of the nonlinear model, as described in (20): The non-principal diagonal elements of its RGA have a value of 1.22, and its stationary condition number is 2.5.Additionally, there is a non-minimum phase multivariable zero with a time constant of 169.78 s.The presence of the multivariable RHP zero imposes control limitations on the achievable bandwidth.To qualitatively validate the model accuracy, an open-loop experiment was conducted around the operating point using the experimental plant.To achieve this, pulses with an amplitude of 10 cm 3 /s are applied to the system inputs, starting from the operating point.Figure 7 compares the real data with the simulation performed using both linear and nonlinear models.Qualitatively, a good fit is observed.To quantitatively measure the fit, the root mean square error (RMSE) value of h 1 and h 2 is used, as defined in (26).These values are presented in Table 3.The results validate both the nonlinear and linear models.

RMSE(h
Actuators 2024, 13, x FOR PEER REVIEW 14 of 22 ( The expressions obtained in ( 17) and ( 21) from linearizing the nonlinear model developed for a manual valve configuration can be used to directly obtain the parameters of the linear models around other operating points without the need to perform any experiments as would be the case with the black-box approach based on identification.Table 4 lists the linear model parameters in (20) for different operating points.

Proposed Control Systems using Decoupling
This section presents three multivariable control strategies for water levels in the lower tanks, starting from the linear model in (25).The strategies highlight the interaction and multivariable zero RHP problems.The proposed controllers are a multiloop PI control, a scheme that combines a multiloop PI control and a static decoupling network, and a centralized control that combines a multiloop PI control with simplified decoupling.Section 5 will present the results of each design, both in simulation and experimentally.The expressions obtained in (17) and (21) from linearizing the nonlinear model developed for a manual valve configuration can be used to directly obtain the parameters of the linear models around other operating points without the need to perform any experiments as would be the case with the black-box approach based on identification.Table 4 lists the linear model parameters in (20) for different operating points.

Proposed Control Systems Using Decoupling
This section presents three multivariable control strategies for water levels in the lower tanks, starting from the linear model in (25).The strategies highlight the interaction and multivariable zero RHP problems.The proposed controllers are a multiloop PI control, a scheme that combines a multiloop PI control and a static decoupling network, and a centralized control that combines a multiloop PI control with simplified decoupling.Section 5 will present the results of each design, both in simulation and experimentally.

Multi-Loop PI Controller
The proposed control system in this section follows the decentralized control scheme depicted in Figure 1a.Since the non-diagonal elements of the process RGA are equal to 1.2, the water levels h 1 and h 2 , which are controlled variables, are cross-paired with the flow rates f R and f L , respectively, which are the manipulated variables.The PI controllers in both loops are then tuned using a phase margin of 60 • as the specification.PI controllers with the parallel structure given in (27) are used, where K P is the proportional gain (cm 2 /s), K I is the integral gain (cm 2 /s 2 ), and h k_ref are level references.When designing, it is important to consider that loop controllers may interact with each other.In this work, we used the iterative methodology proposed in [8], which adjusts the PI controllers in each iteration using single-input, single-output (SISO) methodologies for the equivalent open-loop process (EOP).This methodology considers the interaction with the other loop and the controller calculated in the previous step.The PI parameters for the h 2 -f L loop are K P2 = 3.96 and K I2 = 9.426 × 10 −3 .For the h 1 -f R loop, the PI parameters are K P1 = 4.92 and K I1 = 11.39 × 10 −3 .

Static Decoupling
This control approach follows the control scheme shown in Figure 1c and involves designing a static decoupling network from the process' stationary state gain matrix G lin (0) to reduce interaction at low frequencies [34].The decoupling network is given using (28), which is the inverse of G lin (0).Two PI controllers are then designed using the same methodology as above for the resulting apparent processes from the process and the static decoupler, i.e., the matrix Q(s) = G lin (s)D 0 .The multi-loop PI control employs diagonal pairing and a phase margin of 60 • as specifications for both loops.The resulting PI parameters are K P1 = 1.02 and K I1 = 2.88 × 10 −3 in the first loop, and K P2 = 0.99 and K I2 = 2.73 × 10 −3 in the second loop.

Simplified Decoupling
This decoupling control approach uses a dynamic decoupler to minimize interaction across all frequencies.The simplified decoupling method is employed, as described in (2), and is represented by the transfer matrix D S (s) in (29).The resulting apparent process Q(s) is a diagonal matrix, as shown in (30), and does not exhibit any interaction.Two PI controllers are independently tuned for the corresponding elements of the resulting diagonal apparent process Q(s) = G lin (s)•D S (s), where the multivariable RHP zero of the process appears.The SISO methodology described in [35] is used with a phase margin of 60 • as a specification in both loops.The resulting PI parameters are K P1 = −1.09and K I1 = −1.873× 10 −3 in the first loop, and K P2 = −1.23 and K I2 = −2.15× 10 −3 in the second loop.

Results
This section applies the three designed control systems for both the nonlinear model and the real system.The simulated results are compared with the experimental results for each design and among the three designs.The test carried out in all three cases involves starting from levels of 18 and 19 cm in h 1 and h 2 , respectively, and performing a step change of 4 cm at 80 s in the h 1 reference and at 4600 s in the h 2 reference.The test ends at 8000 s.
Figures 8-10 display the simulation and experimental results for the designed control systems: decentralized, static decoupling, and simplified decoupling, respectively.Each figure presents the main variables of the control system on the left, including the levels of h 1 and h 2 , and the total flow rates supplied by each pump, f L and f R , from top to bottom.On the right, secondary variables such as the tank levels h 3 and h 4 , and the operation percentage of the pumps, p L and p R , are shown.Qualitatively, a good fit between the simulated and experimental data is observed in all three cases, particularly in the primary variables on the left.However, an offset between the simulated and experimental data is observed in the secondary variables on the right.The imprecise adjustment of the multiturn outlet valves may be the cause of the issues with the levels h 3 and h 4 , whereas nonlinearities in the pumps may be responsible for the offset in the pump signals.To address this issue in the pumps, secondary loops were implemented in the QTS to regulate the flow rate of each pump based on the control signals of the proposed control systems.Thus, the absolute values of the percentage of pump usage are not crucial in the proposed model.Table 5 displays the RMSE values of the system variables between the simulation and experimental results.These values quantitatively assess the degree of fit achieved with the proposed nonlinear model.The values are low enough to validate the proposed nonlinear model.The experimental results were used to compare the three proposed control schemes.Table 6 displays the integrated absolute error (IAE) of tank levels 1 and 2 in relation to their references.The total value is calculated using all experimental data, while the tracking value is calculated using only the values corresponding to reference changes in that water level, providing information on the IAE of reference tracking.The decoupling value is calculated using data corresponding to changes in the reference of the other loop, providing information on the decoupling performance.Decentralized control achieves a fast response to changes in the reference with a settling time of approximately 2000 s.However, it exhibits a significant interaction in the opposite loop when such changes occur.Therefore, testing the use of decoupled control schemes was justified.Static decoupling control shows a similar reference tracking response, albeit slightly slower, as reflected by the higher IAE values for tracking.However, the decentralized control fails to achieve the same level of decoupling as loop 2, which has an IAE value that is 35% lower.Furthermore, loop 2 has lower TV values, particularly in fR, which is 28% lower.Static decoupling cannot eliminate all interactions because it is not a dynamic decoupler.The simplified decoupling control achieves nearly perfect decoupling by eliminating interaction in both loops.Its small IAE values for decoupling are 86% lower than those of the other control schemes in loop 1, and 76% and 65% lower in loop 2 compared to decentralized control and static decoupling, respectively.The decoupling process results in a slower response time, with a settling time of approximately 3000 s.This is evidenced by the higher IAE values for tracking, which are approximately twice as high as those of decentralized control and 60% higher than those of control with static decoupling.The inferior response to reference tracking is a result of the constraint imposed by the multivariable RHP zero, which appears in the apparent process elements Q(s) when simplified decoupling is applied.This RHP zero is responsible for the inverse response to changes in the reference.In this QTS configuration, there is a trade-off between the peak of the inverse response and the interaction in the other loop due to the RHP zero.Regarding TV values, the simplified decoupling method has the lowest values, which are approximately 70% lower than those of the decentralized control.The experimental results were used to compare the three proposed control schemes.Table 6 displays the integrated absolute error (IAE) of tank levels 1 and 2 in relation to their references.The total value is calculated using all experimental data, while the tracking value is calculated using only the values corresponding to reference changes in that water level, providing information on the IAE of reference tracking.The decoupling value is calculated using data corresponding to changes in the reference of the other loop, providing information on the decoupling performance.Table 6 displays the total rate of variation (TV) values of the control signals, which are related to the control effort.The expressions for IAE and TV are presented in (31).Decentralized control achieves a fast response to changes in the reference with a settling time of approximately 2000 s.However, it exhibits a significant interaction in the opposite loop when such changes occur.Therefore, testing the use of decoupled control schemes was justified.Static decoupling control shows a similar reference tracking response, albeit slightly slower, as reflected by the higher IAE values for tracking.However, the decentralized control fails to achieve the same level of decoupling as loop 2, which has an IAE value that is 35% lower.Furthermore, loop 2 has lower TV values, particularly in f R , which is 28% lower.Static decoupling cannot eliminate all interactions because it is not a dynamic decoupler.The simplified decoupling control achieves nearly perfect decoupling by eliminating interaction in both loops.Its small IAE values for decoupling are 86% lower than those of the other control schemes in loop 1, and 76% and 65% lower in loop 2 compared to decentralized control and static decoupling, respectively.The decoupling process results in a slower response time, with a settling time of approximately 3000 s.This is evidenced by the higher IAE values for tracking, which are approximately twice as high as those of decentralized control and 60% higher than those of control with static decoupling.The inferior response to reference tracking is a result of the constraint imposed by the multivariable RHP zero, which appears in the apparent process elements Q(s) when simplified decoupling is applied.This RHP zero is responsible for the inverse response to changes in the reference.In this QTS configuration, there is a trade-off between the peak of the inverse response and the interaction in the other loop due to the RHP zero.
Regarding TV values, the simplified decoupling method has the lowest values, which are approximately 70% lower than those of the decentralized control.

Conclusions
This paper presents a practical grey-box modeling procedure for the quadruple-tank system used at the University of Cordoba.The proposed method overcomes the unrealistic assumptions made in other models, such as assuming linear behavior of the actuators and three-way valves.This results in a more realistic, configurable, and reproducible model.The proposed grey-box approach highlights the following points, which can also be useful for modeling other systems with liquid tanks:

•
Secondary control loops are suggested to regulate the flow rate delivered by each pump, making the new control variables the references for those flow rates.The controllers have integral action, resulting in zero stationary error.In addition, the response of these loops is much faster than the dynamics of the tanks, allowing for the requested changes to the flow rate to be almost instantaneous.This approach helps to avoid pump nonlinearities.• The proposed method for modeling the outflows of the tanks is based on level data measured during the experience of emptying the tank and is expressed using the Bernoulli equation.• The percentage of flow diverted to the lower tanks from the three-way valves is modeled using planes as a function of inlet flow and valve position, based on tank filling experiences.To facilitate modeling of these valves, it is useful for the outlet pipes of the three-way valve to have the same length in the construction of the real plant in order to equalize the head losses of both outlet sections.Thus, the distribution of flow is less affected by the total inlet flow, depending on the valve position.
The proposed grey-box approach has resulted in the development of a nonlinear model and has been developed so that it can be adjusted based on valve positions to exhibit little or large interactions, or RHP zero.This model has allowed for the derivation of general expressions of the linearized process model around a generic operating point.These expressions provide an accurate linear model for designing controls at different operating points without the need for new experiments.In contrast, the black-box approach requires experiments at each operating point to obtain a linear model.In addition, the generic expressions of the linear model of the process allow easy implementation of more advanced control strategies such as adaptive control using gain scheduling.For this study, a configuration with a large interaction and RHP zero was selected.Different control strategies were designed based on the linearized model at one operating point and applied to both the nonlinear model and the real process.The validity of the nonlinear model is demonstrated by comparing the simulations and experimental results.The results and conclusions obtained from applying different control designs in the simulation are extrapolated to the real system.
Therefore, it can be concluded that both the real system and the nonlinear model are valuable tools for research and teaching.They offer the opportunity to simulate and experiment with highly interacting processes, including their control systems.The system interactions provide opportunities to test new control theories for robustness and performance with reasonable time and effort.This study presents three control schemes: a decentralized PI control, a static decoupled control, and a simplified decoupled control.We investigated the bandwidth limitation of the control system caused by a multivariable RHP zero when applying perfect decoupling, such as the simplified decoupled control, in both loops.The reference tracking response of the decentralized control is faster, but it comes at the cost of showing large interaction in the opposite loop.It is worth noting that this plant has been successfully used for practical work by students of control courses, as well as by our research group.

Figure 3 .
Figure 3. (a) Schematic of the quadruple-tank system in the laboratory; (b) real quadruple-tank process in the laboratory.

Figure 3 .
Figure 3. (a) Schematic of the quadruple-tank system in the laboratory; (b) real quadruple-tank process in the laboratory.

Figure 4 .
Figure 4. Experimental measurements of the flow fraction to the lower tanks as a function of the total flow rate of each branch and the position of the three-way valve: (a) right branch and (b) left branch.

Figure 4 .
Figure 4. Experimental measurements of the flow fraction to the lower tanks as a function of the total flow rate of each branch and the position of the three-way valve: (a) right branch and (b) left branch.

Figure 5 .
Figure 5. Real emptying curves and their adjustment for valves v1, v2, v3, and v4 when the degree of opening is 3 turns.

Figure 5 .
Figure 5. Real emptying curves and their adjustment for valves v 1 , v 2 , v 3 , and v 4 when the degree of opening is 3 turns.

Figure 6 .
Figure 6.Values of α and β of the valves v1, v2, v3, and v4 for different opening turns and their polynomial fitting.

Figure 6 .
Figure 6.Values of α and β of the valves v 1 , v 2 , v 3 , and v 4 for different opening turns and their polynomial fitting.

Figure 7 .
Figure 7. Open-loop comparison of the linear and nonlinear models with the real process.

Figure 7 .
Figure 7. Open-loop comparison of the linear and nonlinear models with the real process.

Figure 8 .
Figure 8. Experimental and simulated decentralized control results.

Figure 9 .
Figure 9. Experimental and simulated static decoupling control results.

Figure 10 .
Figure 10.Experimental and simulated simplified decoupling control results.

Table 1 .
Proportion γ L (%) of the flow rate f L leading to tank 1 for different values of f L and different positions of V L .

Table 2 .
Proportion γ R (%) of the flow rate f R leading to tank 1 for different values of f R and different positions of V R .

Table 1 .
Proportion γL (%) of the flow rate fL leading to tank 1 for different values of fL and different positions of VL.

Table 2 .
Proportion γR (%) of the flow rate fR leading to tank 1 for different values of fR and different positions of VR.

Table 3 .
RMSE values for the system level variables in the open-loop comparison.

Table 4 .
Parameters of the linearized model for different operation points.

Table 3 .
RMSE values for the system level variables in the open-loop comparison.

Table 4 .
Parameters of the linearized model for different operation points.

Table 6
displays the total rate of

Table 6 .
Performance indices for the different proposed control systems.

Table 5 .
RMSE values for the variables of the system in the different closed-loop systems.

Table 6 .
Performance indices for the different proposed control systems.