Design of Multivariable PID Control Using Iterative Linear Programming and Decoupling

: The design of multivariable process control systems is specially complicated when there are strong interactions between the different control loops, and even more with multiple time delays. This paper proposes an iterative design method of centralized proportional-integral-derivative (PID) controllers for stable linear systems. The methodology is based on the linear parameterization of equivalent loop transfer functions (ELTFs) for centralized control. These functions capture the dynamics of the other loops and, from a prior design, allow solving the design problem at each iteration with linear programming that shapes the Nyquist plot of the ELTFs in the frequency domain, which also avoids the need for approximations. Two optimizations are proposed: (I) maximizing integral gains by fulfilling linear robustness margins in each ELTF and (II) maximizing linear robustness margins by fulfilling minimum bandwidths in each loop. In both optimizations, static decoupling and decoupling at a frequency close to the bandwidth of each loop are included as constraints, which improves the decoupling performance and the procedure convergence. The effectiveness of the method is verified in three simulation examples (square and non-square) and a lab experimental process. The proposed designs achieve a similar or better response when compared to that achieved by other authors.


Introduction
For decades, proportional-integral-derivative (PID) controllers have been the most widespread in industrial processes, where it is estimated that they are currently present in 90% of control loops due to their efficiency and simplicity of implementation [1].Although predictive control (MPC) is now more widespread in the industry, its bandwidth limitations mean that in the regulation layer, PID controllers are still preferred, with MPC being implemented in a higher layer [2].PID controller tuning methods are mainly developed for single-variable loops [3,4].However, most industrial processes are multivariable and have multiple inputs and multiple outputs (MIMOs).For the performance of these systems not to deteriorate, the possible interaction between the inputs and outputs of the process or their possible multiple time delays must be considered.These aspects render these systems more difficult to control and render the design of MIMO control systems more complicated than in the single-variable case [5].To deal with multivariable systems, there are two main approaches: multiloop (or decentralized) control and centralized control.
In multiloop control, the MIMO system is decomposed into several single-variable loops where each output is paired with an input and a controller is designed for each loop, usually a PID controller [6,7].The resulting control is a diagonal matrix controller.
There are different methods for its design, such as detuning, the sequential loop closing method [8], optimization [9,10], direct design [6], or an iterative method based on the concept of effective open loop processes and robustness specifications in the frequency domain [11,12].In cases with less interaction, even PID SISO design methods [13,14] can be directly used.These approaches work well for weakly coupled processes but may be inadequate for highly interacting and diagonally nondominant processes.In these cases, the centralized control approach is more advisable, which can be realized either by combining a decoupler with a diagonal controller [15] or by designing a full matrix controller from the above blocks or directly [16,17].
In the purely centralized scheme, the full matrix controller is the only block responsible for reducing the interactions and controlling the different outputs of the system.The design of these controls can be approached with non-PID elements in the controller matrix and using different methods such as decoupling control [18], optimization procedures [19], H ∝ design [20], and so on.The resulting controllers can be highly complex and even irrational; therefore, several authors have subsequently proposed to approximate the elements of these controls to the PID structure [21], thus obtaining multivariable PID control.Other methodologies have been developed to directly obtain a complete matrix of PID controllers.
There is a large body of work on multivariable PID control design [22] in both the time and frequency domain.Among the design methods in the time domain, iterative methods based on linear matrix inequalities (LMIs) stand out [23][24][25].However, in the case of systems with multiple time delays, the design of PID controllers can become difficult; therefore, frequency domain-based techniques are often used [2].Several methods use a matrix representation of transfer functions and attempt to obtain an approximation of equivalent transfer functions (ETFs) to design multivariable PID controls [26,27].In [28], the approximate ETFs are derived from the effective relative gain array (RGA) of the process.In [29], ETFs were obtained from the relative normalized gain array (RNGA), which approximates the inverse of the process matrix.In [30], they are calculated from the RGA and the average residence time array.Most centralized control methods are designed for square systems, where the number of process inputs and outputs is the same.For non-square systems, the design can be complicated by the uniqueness of the transfer function matrix of the process.One option is to create the process square by adding or removing the input or output variables.Adding outputs impacts the cost of measurement instrumentation, and removing inputs changes the properties of the original system [31].In [32], the ETF concept was extended to non-square systems and proposed multivariable PID control with optimization via particle swarm optimization.In [33], the design method based on ETFs and the RNGA was extended to non-square systems.
Methods based on transfer functions or ETFs tend to become conservative when the process size increases because they require different approximations.To avoid these approximations, some authors have used internal model control [34] or proposed an LMIbased optimization [35].On the other hand, other authors have avoided approximations during their design by directly working with the matrix representation of the process frequency response and using frequency domain specifications.In [36], the controller was designed using the process response at low frequencies to approximate the closed-loop response to that of a reference model.The method was presented for non-square systems based on the pseudoinverse evaluation of the process transfer function matrix at a lowfrequency point.In [37], the frequency response was used to retune a PID MIMO controller.Other works [2,38] proposed the design of centralized PID control based on the frequency response of the process as a convex optimization problem with LMIs.
In our previous paper [39], we proposed an iterative method to design multivariable PID controls for stable and square linear systems based on their frequency response and the concept of equivalent loop transfer function (ELTF).In that study, ELTFs enabled the calculation of the elements for each column of the controller matrix separately using a linear optimization problem if the linear margin proposed in [40] was used as a robustness specification.Such a margin ensures the minimum values of other well-known margins, such as gain margin, phase margin, or maximum sensitivity.Moreover, such a method does not require any approximation when using the frequency response representation.In this paper, we present a new iterative multivariable PID control design methodology based on the previous one, which is extended.Now, instead of designing each column of the controller matrix separately, the entire design is performed in each iteration in a single linear programming problem.For this purpose, the ELTFs are expressed as a function of all the elements of the controller matrix.Thus, a better convergence is sought.In addition to static decoupling, which reduces the interaction at low frequencies, we add the option of including a dynamic decoupling constraint in each loop at a specific frequency close to the desired bandwidth, which helps to render the resulting system diagonally nondominant without imposing such strong constraints as full dynamic decoupling.In the previous study, the optimization problem maximized the integral gains subject to given robustness values; now, the optimization is also posed to maximize the linear robustness margins of each loop while ensuring the given bandwidths.The effectiveness of the method is verified via a simulation of two square linear processes with time delays with dimensions of 2 × 2 and 3 × 3 and a non-square 2 × 3 system.The proposed methods are also tested on a real quadruple tank process.Comparison with other authors shows that the proposed methodologies achieve a similar or better performance.
The rest of the article is as follows: Section 2 presents the linear robustness margin and the concept of an ETLF in centralized control.Section 3 describes the proposed methodology, explaining the different types of linear programming according to the cost function to be maximized and the different constraints to which it may be subject.In Section 4, the proposed method is applied to different processes and compared with controllers developed by other authors.Finally, Section 5 summarizes the main conclusions.

Preliminaries
Figure 1 shows a centralized unitary feedback control system in which the linear multivariable process is represented by the transfer function matrix G(s) of the dimension n × m and the centralized control by the matrix K(s) of m × n dimensions.When the number of outputs and inputs of the process is the same, these matrices are square.Each element g ij (s) of G(s) represents the transfer function between output y i and input u j .Each controller k ij (s) represents the control element between the control signal u i and the error signal e j , which is the difference between the reference signal r j and the controlled output y j .The open-loop transfer matrix is denoted as L(s) = G(s)K(s) when unitary feedback is assumed.In the 2 × 2 case, the elements l ij (s) are provided in (1), where the Laplace variable s has been omitted, as shown in Figure 1.The elements of a single column of L(s) depend on the elements of a single column of K(s).Something similar happens for higher dimensional processes.Thus, for example, the generic element l ij (s) of the j-th column of L(s) for n × n processes are provided with (2) and depends on the elements of the j-th column of K(s).
the concept of equivalent loop transfer function (ELTF).In that study, ELTFs enabled the calculation of the elements for each column of the controller matrix separately using a linear optimization problem if the linear margin proposed in [40] was used as a robustness specification.Such a margin ensures the minimum values of other well-known margins, such as gain margin, phase margin, or maximum sensitivity.Moreover, such a method does not require any approximation when using the frequency response representation.In this paper, we present a new iterative multivariable PID control design methodology based on the previous one, which is extended.Now, instead of designing each column of the controller matrix separately, the entire design is performed in each iteration in a single linear programming problem.For this purpose, the ELTFs are expressed as a function of all the elements of the controller matrix.Thus, a better convergence is sought.In addition to static decoupling, which reduces the interaction at low frequencies, we add the option of including a dynamic decoupling constraint in each loop at a specific frequency close to the desired bandwidth, which helps to render the resulting system diagonally nondominant without imposing such strong constraints as full dynamic decoupling.In the previous study, the optimization problem maximized the integral gains subject to given robustness values; now, the optimization is also posed to maximize the linear robustness margins of each loop while ensuring the given bandwidths.The effectiveness of the method is verified via a simulation of two square linear processes with time delays with dimensions of 2 × 2 and 3 × 3 and a non-square 2 × 3 system.The proposed methods are also tested on a real quadruple tank process.Comparison with other authors shows that the proposed methodologies achieve a similar or better performance.
The rest of the article is as follows: Section 2 presents the linear robustness margin and the concept of an ETLF in centralized control.Section 3 describes the proposed methodology, explaining the different types of linear programming according to the cost function to be maximized and the different constraints to which it may be subject.In Section 4, the proposed method is applied to different processes and compared with controllers developed by other authors.Finally, Section 5 summarizes the main conclusions.

Preliminaries
Figure 1 shows a centralized unitary feedback control system in which the linear multivariable process is represented by the transfer function matrix G(s) of the dimension n × m and the centralized control by the matrix K(s) of m × n dimensions.When the number of outputs and inputs of the process is the same, these matrices are square.Each element gij(s) of G(s) represents the transfer function between output yi and input uj.Each controller kij(s) represents the control element between the control signal ui and the error signal ej, which is the difference between the reference signal rj and the controlled output yj.The open-loop transfer matrix is denoted as L(s) = G(s)K(s) when unitary feedback is assumed.In the 2 × 2 case, the elements lij(s) are provided in (1), where the Laplace variable s has been omitted, as shown in Figure 1.The elements of a single column of L(s) depend on the elements of a single column of K(s).Something similar happens for higher dimensional processes.Thus, for example, the generic element lij(s) of the j-th column of L(s) for n × n processes are provided with (2) and depends on the elements of the j-th column of K(s).
The problem with these open-loop transfer functions is that they do not capture the interactions and dynamics with the other loops when they are closed.In [39], the concept of an equivalent open-loop transfer function (ELTF) for centralized control was introduced.This function represents the transfer function between an output and its error signal when the loop is open and the rest of the loops remain closed.Figure 2 shows the ELTF l 1 (s) framed for the 2 × 2 case.Unlike the l ij (s) functions in L(s), the ELTF incorporates the interactions and dynamics with the other control loops.For the case with two outputs, the expressions of the ELTFs are shown in (3).For n dimensions, the generic expression is provided in (4).
Electronics 2024, 13, x FOR PEER REVIEW 4 of 27 The problem with these open-loop transfer functions is that they do not capture the interactions and dynamics with the other loops when they are closed.In [39], the concept of an equivalent open-loop transfer function (ELTF) for centralized control was introduced.This function represents the transfer function between an output and its error signal when the loop is open and the rest of the loops remain closed.Figure 2 shows the ELTF l1(s) framed for the 2 × 2 case.Unlike the lij(s) functions in L(s), the ELTF incorporates the interactions and dynamics with the other control loops.For the case with two outputs, the expressions of the ELTFs are shown in (3).For n dimensions, the generic expression is provided in (4). ) Figure 2. Equivalent loop transfer function l1(s) for 2 × 2 processes [39].
As demonstrated in [39], ELTFs can be expressed as a function of controllers in the same column, which can be designed in isolation in an iterative process if the other control columns are known.By using PID controllers with parallel structures and linear robustness margins as a design specification for the frequency response of the ELTFs, this problem was posed as linear programming for each ELTF.The linear margin Lm is shown in Figure 3, where it is determined via the straight line r intercepting the negative real axis at the point (−1 + Lm, 0) with an angle α.If the ELTF li(s) is stable, closed-loop stability is guaranteed if its Nyquist trace is below the r-straight line because it does not encircle the critical point −1.Using the frequency response li(jω), where ω is the frequency variable, this condition can be expressed with the inequation shown in (5), where ℑ() represents the imaginary part, ℜ() is the real part, and tan() is the tangent function.Furthermore, as the polar trace lies below the line r, the minimum values of other classical robustness margins, such as gain margin Am, phase margin ϕm, or maximum sensitivity margin Ms, all depicted in Figure 3, are guaranteed.More details on the relationship between these classical robustness margins and the values of Lm and α can be found in [39,40].The robustness As demonstrated in [39], ELTFs can be expressed as a function of controllers in the same column, which can be designed in isolation in an iterative process if the other control columns are known.By using PID controllers with parallel structures and linear robustness margins as a design specification for the frequency response of the ELTFs, this problem was posed as linear programming for each ELTF.The linear margin L m is shown in Figure 3, where it is determined via the straight line r intercepting the negative real axis at the point (−1 + L m , 0) with an angle α.If the ELTF l i (s) is stable, closed-loop stability is guaranteed if its Nyquist trace is below the r-straight line because it does not encircle the critical point −1.Using the frequency response l i (jω), where ω is the frequency variable, this condition can be expressed with the inequation shown in (5), where ℑ() represents the imaginary part, ℜ() is the real part, and tan() is the tangent function.Furthermore, as the polar trace lies below the line r, the minimum values of other classical robustness margins, such as gain margin A m , phase margin ϕ m , or maximum sensitivity margin M s , all depicted in Figure 3, are guaranteed.More details on the relationship between these classical robustness margins and the values of L m and α can be found in [39,40].The robustness increases with increasing values of L m and α, with typical ranges of [0.5, 0.85] for L m and [35 • , 90 • ] for α.Moreover, the higher these values are, the more damped the system will be and the lower the control effort will be; however, as a trade-off, the bandwidth decreases, which tends to render the settling time of the time response longer.

ℑ(l
increases with increasing values of Lm and α, with typical ranges of [0.5, 0.85] for Lm and [35°, 90°] for α.Moreover, the higher these values are, the more damped the system will be and the lower the control effort will be; however, as a trade-off, the bandwidth decreases, which tends to render the settling time of the time response longer.
Figure 3. Linear margin of robustness in the Nyquist plot [12].
In addition, Figure 3 shows the line q, which is tangential to the unit circle in the third quadrant and intercepts the negative real semi-axis with an angle β.With this line, a condition can be established on li(jω) to ensure a lower limit ωx at the crossover frequency ωc, which marks the bandwidth of the loop.Such a condition is provided in the inequations in (6): the first forces the trace of li(jω) to be below line q for frequencies lower than ωx, and the second forces the trace at higher frequencies to be above line q.Here, cos and sin are the sine and cosine functions, respectively.

Proposed Methodology
This section describes the proposed iterative method for the design of MIMO PID controllers based on the frequency response of the ELTFs and the linear margin.The proposed method performs the design of all PID controllers in each iteration and not separately by columns, as in our previous work [39]; moreover, the linear programming developed for its resolution is proposed with different cost functions and constraints.

Linear Parameterization of ELTFs
The proposed method is based on a linear parameterization of the ELTFs as a function of the PID controller parameters.Assuming a parallel PID structure in the controller matrix K(s), the controllers have the form provided in (7), where KPij is the proportional gain of the control element kij(s), KIij is the integral gain, and KDij is the derivative gain.These can be linearly parameterized as a function of the parameter vector ρij with the above gains, as shown in (7).In addition, Figure 3 shows the line q, which is tangential to the unit circle in the third quadrant and intercepts the negative real semi-axis with an angle β.With this line, a condition can be established on l i (jω) to ensure a lower limit ω x at the crossover frequency ω c , which marks the bandwidth of the loop.Such a condition is provided in the inequations in (6): the first forces the trace of l i (jω) to be below line q for frequencies lower than ω x , and the second forces the trace at higher frequencies to be above line q.Here, cos and sin are the sine and cosine functions, respectively.

Proposed Methodology
This section describes the proposed iterative method for the design of MIMO PID controllers based on the frequency response of the ELTFs and the linear margin.The proposed method performs the design of all PID controllers in each iteration and not separately by columns, as in our previous work [39]; moreover, the linear programming developed for its resolution is proposed with different cost functions and constraints.

Linear Parameterization of ELTFs
The proposed method is based on a linear parameterization of the ELTFs as a function of the PID controller parameters.Assuming a parallel PID structure in the controller matrix K(s), the controllers have the form provided in (7), where K Pij is the proportional gain of the control element kij(s), K Iij is the integral gain, and K Dij is the derivative gain.These can be linearly parameterized as a function of the parameter vector ρ ij with the above gains, as shown in (7).
We defined a vector ρ with all the parameters of all the K(s) controllers ordered by rows as shown in (8) for the 2 × 2 case.
We defined the product of a process element g km (s) by a controller element k ij (s), as shown in (9), where the vector φ km (s) encompasses the process element g km (s) and the vector [1, 1/s, s] associated with the PID element.This product was linearly parameterized with respect to the PID parameter vector ρ ij and the three-element row vector φ km (s).For the elements of this vector, a priori, to be positive elements, the corresponding sign was extracted.The signs of the controller elements were obtained, assuming that the centralized control matrix K(s) at low frequencies would be proportional to the inverse (or pseudoinverse) of the steady-state process matrix, i.e., G(0).In fact, the sign was extracted so that the integral gains K Iij were positive since they were forced to be positive in the subsequent optimization.After optimization, the sign corresponding to each controller k ij (s) was returned.The sign function is represented by sign().
From the latter expression, expressed the elements of the matrix L(s) as a linear parameterization of the total PID parameter vector ρ with the row vectors Λ ij (s), which captured the process dynamics corresponding to l ij (s) and the part in s of the PID controllers.Thus, for the 2 × 2 case, provided with the vector ρ in (8), the elements of L(s) in ( 1) could be compactly expressed as follows in (10), where Λ ij (s) and ρ are 12-element vectors.For 3 × 3 or higher cases, a similar procedure could be followed to finally have any element of L(s) linearly parameterized with respect to the vector ρ of the PID parameters as l ij (s) = Λ ij (s)ρ.
From these parameterizations of the open-loop elements l ij (s), the ELTFs could be expressed as a linear parameterization of the PID parameter vector, assuming that some functions of L(s) are previously known.For example, given the generic ELTF expression l j (s) provided in (4), in the summation terms, the denominator 1 + l ii (s) must be assumed to be previously known from a previous design, as well as one of the two elements of the numerator product.For instance, with the expressions in (3) for the 2 × 2 case, the ELTF l 1 (s) could be parameterized as shown in (11), where the functions l 22 (s) and l 12 (s) are assumed to be known and where Ψ 1a (s) is a 12-element row vector.In this case, the columns of Ψ 1a (s) corresponding to the parameters of k 12 and k 22 are zeros, then l 1 (s) is parameterized as a function of only the parameters of the controllers in column 1 of K(s), which would be the vectors ρ 11 and ρ 21 .However, l 1 (s) could also be expressed as in (12) if, instead of assuming l 12 (s) to be known, l 21 (s) is assumed to be known.In this case, l 1 (s) depends on all controllers of K(s), i.e., on all the elements of ρ.For this case, we used the nomenclature Ψ 1b (s).For l 2 (s), a similar parameterization could also be performed with two options.Using these two ways of representing the ELTFs allows better convergence in the subsequently proposed iterative process.
Proceeding similarly for the higher-dimensional processes, the ELTFs could also be expressed as a linear parameterization of the PID parameter vector.Thus, for example, for the 3 × 3 case and from (4), the ELTF l 1 (s) could be expressed as a function of only the controllers in the first column of K(s) with the expression in (13), where columns 2 and 3 of L(s) are assumed to be known.The function l 1 (s) could also be expressed as a function of all the controllers with the expression in (14), where other elements of L(s) are assumed to be known.In these expressions, the PID parameter vector would be a column vector of size 27, and the vectors Λ ij (s), Ψ 1a (s), and Ψ 1b (s) would be row vectors of size 27.In this study, the two linear parametrizations to express an ELTF are (a)

PID Parameter Tuning via Linear Programming
The expressions of the ELTFs seen in the previous section can cause very complex and even irrational functions in the case of processes with time delays.In some design methods, this would imply approximations with the corresponding loss of information.To avoid this, we propose to directly work with the array of the frequency response of the elements of G(s), L(s) and the ELTFs.In this way, each point of the Nyquist trace of the ELTF l i (jω k ) can be expressed as a linear function of the PID parameter vector, as seen in the previous section, and decomposed into real and imaginary parts, as shown in ( 15) for option a.
From these expressions of the ELTFs in the frequency domain using the linear robustness margin and starting from an initial L(s) matrix, the design of the MIMO PID control is proposed as a linear programming problem that determines the vector ρ with all the PID parameters.For this purpose, the objective function and the constraints must be linearly parameterized with respect to the vector ρ.
Two types of optimizations are proposed in this study: (I) an optimization that maximizes the performance of the loops subjected to certain minimum robustness values and (II) another optimization that maximizes the robustness of the loops by ensuring minimum bandwidth values in the loops.Although these optimizations are detailed below for the case of 2 × 2 processes, they are easily scalable to larger systems.Furthermore, although these optimizations are presented to obtain PID controllers in the K(s) matrix, one can easily force some elements to have a PI structure if the corresponding derivative gains K D are constrained to be equal to zero.

Optimization I
This optimization aims to maximize the integral gains K I of the PID controllers in K(s), with the ELTFs being subject to minimum robustness requirements provided with the linear margin defined by the values of L m and α.This is intended to optimize the closed-loop performance in terms of disturbance rejection because the integrated error (IE) of the transient response is minimized when the gains K I are maximized [41].
If the elements of the transfer matrix of the process G(s) are stable and the Nyquist plots of the ELTFs lie below the straight line representing the linear margin, a condition provided in the inequation (5), the closed-loop system is stable.Furthermore, since denominators with the form 1 + l ii (s) appear in the ELTF expressions, as shown in (3) and (4), it is also necessary that such terms have no zeros in the right half-plane; otherwise, they would produce unstable poles in the ELTFs, which would render the linear margin constraint inapplicable.Therefore, it is also necessary to consider as a constraint that the Nyquist plots of the diagonal open-loop elements l ii (s) of L(s) do not surround the critical point −1.Furthermore, the linear margin constraint provided in (5) is applied using the two forms (a and b) of each ELTF to improve the convergence of the method.Substituting (15) into (5), the constraint for the ELTF of loop i and option a is as follows in (16), where cot is the cotangent function: To ensure that the open-loop function l ii (s) does not surround the point −1, we propose to use a constraint similar to the previous one but applied to a parallel line (with the same value of α) but intercepting at −0.8 the negative real semi-axis, which is equivalent to L m = 0.2.A lower value of L m , even L m = 0, could also be used and then, the line would intersect closer to point −1; however, we consider 0.2 to be a good value to increase the robustness to modeling or computational errors in this constraint without imposing a severe limitation on the feasible search space.Such a constraint is linearly parameterized with respect to ρ as follows: Another constraint that enhances convergence and improves decoupling is the imposed static decoupling, i.e., for ω = 0.This implies that the off-diagonal elements of L(j0) should be equal to zero.This constraint is applied only to the factors multiplied by the integral gains as shown in (18) for the 2 × 2 case, since these terms are the most important at low frequencies.This constraint is very useful in systems with weak diagonal dominance.
Another constraint that helps to further improve the decoupling performance and renders the resulting system diagonally dominant is to impose decoupling at a frequency ω x close to the cutoff frequency of the bandwidth.This constraint is best used to try to improve the performance of a previous design created without this constraint and to use the crossover frequencies calculated with that design as the ω x frequencies at which to decouple.Decoupling at frequency ω xj , where j is the control loop, is provided in the following condition in (19), which indicates that at that frequency, the elements l ij (jω xj ) of column j must be zero.This constraint also facilitates faster convergence as the search space is reduced.
The proposed optimization I is expressed as the linear programming problem provided in (20), where the objective function and the constraints are parameterized with respect to the PID parameter vector ρ.The objective function maximizes the integral gains, which are multiplied by unity; the remaining parameters are multiplied by zero.The constraints incorporated in (20) are as follows: (a) Inequation for fulfilling the linear margin by each ELTF to ensure minimum robustness values.To improve the convergence of the iterative method, the two proposed expressions of the ELTFs are used: one that depends only on the controllers of the associated column (option a) and another that depends on all the control elements (option b); (b) Inequation to prevent the Nyquist polar trace of the diagonal elements of the openloop matrix L(s) from encircling the critical point.Thus, the ELTFs are stable if the elements of G(s) are stable; (c) Equality to achieve static decoupling.As an example, the expressions for the 2 × 2 case in (18) are parameterized with respect to ρ.Here, 0 1×3 represents a row vector of three zero elements; (d) Equality to achieve decoupling at the selected ω x frequencies.The constraint in (19) is linearly parameterized with respect to ρ; (e) Lower (LB) and upper (UB) bounds on the design parameters.Integral gains are forced to be positive, whereas proportional and derivative gains can have any sign.

Optimization II
The second proposed optimization attempts to maximize the linear robustness margin L m of each loop with the ELTFs subject to minimum bandwidth requirements, which is related to the crossing frequency ω c of the phase margin, i.e., the frequency at which the Nyquist plot crosses the unit circle.The higher the crossover frequency, the shorter the rise time and the faster the disturbance rejection.
The conditions defined in ( 6) are used to define the bandwidth requirements; in them, the frequency ω x is a lower approximation of the crossover frequency ω c .The user must specify the angle β defining the line q in Figure 3 and the desired frequency ω x .This frequency cannot be too high; otherwise, linear programming will not find a solution.For example, in systems with a time delay, ω x must be less than 5/t d , t d being the longest time delay associated with the ELTF.Substituting ( 15) in (6), such constraints for the ELTF of loop i can be expressed as a line parameterization of the PID parameter vector ρ as follows in (21) for the ELTF represented with option a.This indicates that the Nyquist plot of the ELTF must lie below the line q at frequencies below ω x and above it for higher frequencies. [cos Furthermore, in this optimization, it must be considered that the distances L m of the linear margins of robustness to be maximized in each loop must be added as a variable to the vector of decision variables.Thus, as an example, for the 2 × 2 case, the new design vector θ would be defined as shown in (22): first, the vector of the PID parameters ρ and then the two linear margins are symbolized by ℓ 1 and ℓ 2 .On the latter, a lower bound of 0.3 and a maximum of 0.95 are imposed.
In the case of 2 × 2 processes, the linear programming method defining the proposed optimization II is shown in (23), where the objective function and constraints are parameterized with respect to the new vector of the decision variables θ in (22).The objective function maximizes ℓ 1 and ℓ 2 , which are the only variables multiplied by 1.Some constraints in (23) are similar to those in optimization problem I, but others change, are sometimes not advised, or are new.These constraints are as follows: (a) Inequations to meet the bandwidth specification ω x in each ELTF.To improve the convergence of the iterative method, the two proposed expressions of the ELTFs are used: one that only depends on the associated column controllers (option a) and one that depends on all the control elements (option b); (b) Inequation for the fulfillment of the linear margin by each ELTF similar to optimization I but with two differences: the linear margin is now a decision variable; and, in addition, this constraint is now only required for frequencies higher than ω x .For lower frequencies, it is not required to be to the right of line r in Figure 3 since the previous constraint renders the Nyquist plot of the ELTF to be below the line q defined by β and tangent to the unit circle for frequencies below ω x .This constraint also causes the polar trace to move away from the critical point −1; (c) Inequation to prevent the polar trace of the diagonal elements of the open-loop matrix L(s) from encircling the critical point.It is similar to optimization I but is parameterized with respect to θ; (d) Equality to achieve static decoupling.This is similar to optimization I but is parameterized with respect to θ; (e) Equality to achieve decoupling at frequency ω x .This is similar to optimization I but is parameterized with respect to θ.However, it must be said that in optimization II, adding this may reduce the achievable bandwidth; (f) Lower (LB) and upper (UB) bounds on the design parameters.maximize 0 1x12 1 1 θ subject to :

Iterative Algorithm
The optimizations proposed in the previous section are jointly performed for the ELTFs of the multivariable system to calculate the total PID parameter vector ρ.Because in some constraints some open-loop functions l ij (s) are assumed to be known, an initial control matrix is required to calculate the open-loop matrix L(s).When the new control matrix K(s) is obtained after optimization, it is necessary to recalculate the matrix L(s), since it will have changed and recheck whether the design specifications are met.If not, iteration is necessary.Therefore, the proposed design method is iterative, and the input data are the process transfer matrix G(s), the design specifications, and an initial control matrix K 0 (s).As the output, it provides the resultant multivariable PID control.
Figure 4 shows the flowchart of the proposed iterative algorithm.First, the iteration variable m is set to 0. Next, an initial control matrix K 0 (s) is chosen, and the frequency response of the open-loop matrix L 0 (jω) is calculated before entering the iterative loop.Although the initial K 0 (s) matrix may somewhat affect the number of iterations, in the analyzed examples, practically the same final controls were obtained using different guesses.Here, we propose to use the inverse of the static gain matrix of the process, which is G(0), because of its computational simplicity and good convergence speed.In the case of nonsquare processes, the pseudoinverse of G(0) can be used.Then, we enter an iterative loop in which in each iteration m, the following steps are performed: 1.
Calculate, according to Section 3.1, the vectors Ψ ia (jω k ), Ψ ib (jω k ) for the ELTF of each loop and the vectors Λ ij (jω k ) from the frequency response of the open-loop matrix known from the previous iteration L m−1 (jω); 2.
Obtain the PID parameter vector ρ m using the linear programming method proposed in Section 3.2 depending on the design specifications; If the termination condition is fulfilled, the design is accepted and the resulting vector ρ is used to generate the multivariable PID control matrix K(s); otherwise, the iteration variable m is incremented and the process iterates again.The termination condition is satisfied if, during three consecutive iterations, the variation of the decision variables is below a user-defined tolerance.A maximum number of iterations before the process stops in the case of no convergence is also defined.After calculating the control parameter vector ρ, the PID elements of the control matrix K(s) are implemented.When implementing the designed PID controllers, the PID structure provided in (24) is used, where u(s) is the control signal, y(s) is the controlled variable, and e(s) is the error signal.This PID structure avoids excessive spikes in the control signal caused by derivative action when there are sudden changes in the reference signal [41].Derivative action is filtered with a first-order term to render the controllers realizable.A value of N = 20 is chosen so that the effect of the filter can be neglected in the frequency range of interest.

Illustrative Examples
In this section, the effectiveness of the proposed methodology is verified in a simulation by applying it to three stable processes of different dimensions: a 2 × 2 square system, a 3 × 3 square system, and a 2 × 3 non-square system.In each example, the proposed multivariable PID control is compared with those proposed by other authors.

Example 1: A Wood and Berry Distillation Column
This example uses the distillation column process that is defined by the 2 × 2 transfer matrix GWB(s) in ( 25) [42].It is a classical testbench in multivariate systems because it presents significant interactions and contains multiple time delays.The time constants are expressed in minutes.From GWB(s), a frequency response array of 1000 elements logarithmically spaced in the frequency range of [10 −5 , 5] rad/min is computed.To design the multivariable PID control, optimization II is used with bandwidth design specifications of 0.4 and 0.18 rad/min for loops 1 and 2, respectively.These values were chosen to obtain a response speed similar to that of other authors.In addition, values of α = 70° and β = 35° After calculating the control parameter vector ρ, the PID elements of the control matrix K(s) are implemented.When implementing the designed PID controllers, the PID structure provided in (24) is used, where u(s) is the control signal, y(s) is the controlled variable, and e(s) is the error signal.This PID structure avoids excessive spikes in the control signal caused by derivative action when there are sudden changes in the reference signal [41].Derivative action is filtered with a first-order term to render the controllers realizable.A value of N = 20 is chosen so that the effect of the filter can be neglected in the frequency range of interest.

Illustrative Examples
In this section, the effectiveness of the proposed methodology is verified in a simulation by applying it to three stable processes of different dimensions: a 2 × 2 square system, a 3 × 3 square system, and a 2 × 3 non-square system.In each example, the proposed multivariable PID control is compared with those proposed by other authors.

Example 1: A Wood and Berry Distillation Column
This example uses the distillation column process that is defined by the 2 × 2 transfer matrix G WB (s) in (25) [42].It is a classical testbench in multivariate systems because it presents significant interactions and contains multiple time delays.The time constants are expressed in minutes.From G WB (s), a frequency response array of 1000 elements logarithmically spaced in the frequency range of [10 −5 , 5] rad/min is computed.To design the multivariable PID control, optimization II is used with bandwidth design specifications of 0.4 and 0.18 rad/min for loops 1 and 2, respectively.These values were chosen to obtain a response speed similar to that of other authors.In addition, values of α = 70 • and β = 35 • are specified for both loops, and the two decoupling constraints (d and e) shown in (23) are used.By using the inverse of G WB (0) as the initial controller, the proposed method obtains the multivariable PID control shown in (26) in only five iterations, as shown in Figure 5.This figure shows the progression of the PID parameters of each controller k ij (s) and the margins L m during the iterative procedure.Considering that the proposed procedure always involves at least three iterations to ensure convergence, in this example, the convergence is quickly achieved.This is partly due to the use of all the proposed constraints for optimization II described in Section 3.2.2.The linear margin values achieved are ℓ 1 = 0.721 and ℓ 2 = 0.704.
By using the inverse of GWB(0) as the initial controller, the proposed method obtains the multivariable PID control shown in (26) in only five iterations, as shown in Figure 5.This figure shows the progression of the PID parameters of each controller kij(s) and the margins Lm during the iterative procedure.Considering that the proposed procedure always involves at least three iterations to ensure convergence, in this example, the convergence is quickly achieved.This is partly due to the use of all the proposed constraints for optimization II described in Section 3.2.2.The linear margin values achieved are ℓ1 = 0.721 and ℓ2 = 0.704.
In the top part of Figure 6, the resulting Nyquist plots for the two ELTFs are shown along with the r and q lines associated with the design constraints.It should be recalled that, in optimization II, the linear margin constraint only applies at frequencies above ωx.Table 1 lists the other classical robustness margins and crossover frequency ωc achieved with the proposed design with those obtained by other authors.Specifically, these other controllers are the centralized PID control designed via an iterative LMI approach by Feng et al. [43] and the centralized PID control for decoupling designed by Jeng et al. via a databased method in [21].The proposed control achieves the highest gain margins Am, the best sensitivity margins Ms, and the highest bandwidth in loop 1.The Ms index is calculated for the corresponding ELTF according to the expression provided in (27).In the top part of Figure 6, the resulting Nyquist plots for the two ELTFs are shown along with the r and q lines associated with the design constraints.It should be recalled that, in optimization II, the linear margin constraint only applies at frequencies above ω x .Table 1 lists the other classical robustness margins and crossover frequency ω c achieved with the proposed design with those obtained by other authors.Specifically, these other controllers are the centralized PID control designed via an iterative LMI approach by Feng et al. [43] and the centralized PID control for decoupling designed by Jeng et al. via a data-based method in [21].The proposed control achieves the highest gain margins A m , the best sensitivity margins M s , and the highest bandwidth in loop 1.The M s index is calculated for the corresponding ELTF according to the expression provided in (27).The bottom part of Figure 6 shows the Nyquist plots of the diagonal elements of the resulting L(s) matrix with the corresponding Gershgorin circles.These results show the diagonal dominance of the resulting open-loop process and how perfect decoupling has been practically achieved at the crossover frequency since the circles around ωc are almost of zero radius.
Figure 7 shows the closed-loop time response of the proposed PID MIMO control compared with the other two above-mentioned designs.The simulations were performed using a unit step change at t = 1 min in the first reference and at t = 100 min in the second reference.In addition, at t = 200 min, there is a load disturbance consisting of a step change of amplitude 0.5 at both process inputs.From the simulation data, the integrated absolute error (IAE) and the control signal total variation (TV) of each loop are calculated in the time domain.The general expressions are shown in (27).These values are used as performance indices of the time response: the first is related to the error and the second is related to the control effort.These values are also illustrated in Table 1.The bottom part of Figure 6 shows the Nyquist plots of the diagonal elements of the resulting L(s) matrix with the corresponding Gershgorin circles.These results show the diagonal dominance of the resulting open-loop process and how perfect decoupling has been practically achieved at the crossover frequency since the circles around ω c are almost of zero radius.
Figure 7 shows the closed-loop time response of the proposed PID MIMO control compared with the other two above-mentioned designs.The simulations were performed using a unit step change at t = 1 min in the first reference and at t = 100 min in the second reference.In addition, at t = 200 min, there is a load disturbance consisting of a step change of amplitude 0.5 at both process inputs.From the simulation data, the integrated absolute error (IAE) and the control signal total variation (TV) of each loop are calculated in the time domain.The general expressions are shown in (27).These values are used as performance indices of the time response: the first is related to the error and the second is related to the control effort.These values are also illustrated in Table 1.
The proposed MIMO PID control shows a settling time similar to Jeng's control in loop 1 and somewhat higher in loop 2. Feng's control has the highest settling time in both loops.The proposed control shows excellent decoupling performance, especially compared to Feng's control.Jeng's control achieves almost perfect decoupling because its design is based on full dynamic decoupling; however, this means its response to load disturbance is the slowest of them all.The proposed PID control rejects the load disturbances earlier than the others, which causes it to have the lowest IAE values in Table 1.The IAE values in loops 1 and 2 obtained with the proposed control are, respectively, 24.2% and 12.9% lower than Feng's values and 33% and 38.2% lower than those of Jeng's control.Moreover, this response is achieved with TV values similar to or lower than those of the other methods.The TV values of the proposed control in loops 1 and 2 are 13.3% and 26.7% lower than Feng's values, respectively.Compared to Jeng's control, the proposed control achieves better IAE values at the expense of an increase of approximately 30% in the TV values of both loops.The proposed MIMO PID control shows a settling time similar to Jeng's control loop 1 and somewhat higher in loop 2. Feng's control has the highest settling time in bo loops.The proposed control shows excellent decoupling performance, especially com pared to Feng's control.Jeng's control achieves almost perfect decoupling because its d sign is based on full dynamic decoupling; however, this means its response to load di turbance is the slowest of them all.The proposed PID control rejects the load disturbanc earlier than the others, which causes it to have the lowest IAE values in Table 1.The IA values in loops 1 and 2 obtained with the proposed control are, respectively, 24.2% an 12.9% lower than Feng's values and 33% and 38.2% lower than those of Jeng's contro Moreover, this response is achieved with TV values similar to or lower than those of th other methods.The TV values of the proposed control in loops 1 and 2 are 13.3% an 26.7% lower than Feng's values, respectively.Compared to Jeng's control, the propose control achieves better IAE values at the expense of an increase of approximately 30% the TV values of both loops.

Example 2: An Ogunnaike and Ray Distillation Column
The transfer matrix of this process is provided in the 3 × 3 GOR(s) matrix in (28) [44 This is a process with multiple time delays and significant interactions, where the tim constants are expressed in minutes.From this matrix, a frequency response array of 100 elements in the range of [10 −5 , 5] rad/min is obtained.In this example, optimization I applied with robustness specifications provided in Lm, with values of 0.85, 0.8, and 0.6 and α values equal to 85°, 80°, and 65° for loops 1, 2, and 3, respectively.In addition, th decoupling constraint d is applied at frequencies 0.16, 0.18, and 0.8 rad/min, respectivel in each loop.After five iterations, the multivariable PID control provided in ( 29) is o tained.Figure 8 shows the progression of the PID parameters using the inverse of th GOR(0) matrix as the initial control matrix K 0 .Although this second example has a high dimension than the previous example, the proposed iterative procedure quickly co verges.

Example 2: An Ogunnaike and Ray Distillation Column
The transfer matrix of this process is provided in the 3 × 3 G OR (s) matrix in (28) [44].This is a process with multiple time delays and significant interactions, where the time constants are expressed in minutes.From this matrix, a frequency response array of 1000 elements in the range of [10 −5 , 5] rad/min is obtained.In this example, optimization I is applied with robustness specifications provided in L m , with values of 0.85, 0.8, and 0.65, and α values equal to 85 • , 80 • , and 65 • for loops 1, 2, and 3, respectively.In addition, the decoupling constraint d is applied at frequencies 0.16, 0.18, and 0.8 rad/min, respectively, in each loop.After five iterations, the multivariable PID control provided in (29) is obtained.Figure 8 shows the progression of the PID parameters using the inverse of the G OR (0) matrix as the initial control matrix K 0 .Although this second example has a higher dimension than the previous example, the proposed iterative procedure quickly converges.
The top part of Figure 9 shows the Nyquist plots of the three resulting ELTFs, where the fulfillment of the linear margin constraint is verified.Table 2 shows the other classical robustness indices for each loop.For comparison, the indices obtained with the controllers proposed by other authors are also shown.One of these controllers is the decentralized PID control of Silva and Barros [45], which was designed via an iterative method in the frequency domain.The other is the centralized PI control of Ram and Chidambaram [17], which is based on the combination of static decoupling and multiloop PI control.The proposed control achieves very robust margins mainly in the first two loops, where it has the largest phase margins and the lowest maximum sensitivity values; however, it has the lowest cutoff frequencies.Ram's proposed control has the most robust design for the third loop.The bottom part of Figure 9 shows the Nyquist plots of the diagonal open-loop processes of L(s) with the Gershgorin circles.The small circles in l 11 (s) and l 22 (s) indicate good diagonal dominance in these loops, not so in l 33 (s) with very large circles, which causes this loop to exhibit significant interaction.
The top part of Figure 9 shows the Nyquist plots of the three resulting ELTFs, where the fulfillment of the linear margin constraint is verified.Table 2 shows the other classical robustness indices for each loop.For comparison, the indices obtained with the controllers proposed by other authors are also shown.One of these controllers is the decentralized PID control of Silva and Barros [45], which was designed via an iterative method in the frequency domain.The other is the centralized PI control of Ram and Chidambaram [17], which is based on the combination of static decoupling and multiloop PI control.The proposed control achieves very robust margins mainly in the first two loops, where it has the largest phase margins and the lowest maximum sensitivity values; however, it has the lowest cutoff frequencies.Ram's proposed control has the most robust design for the third loop.The bottom part of Figure 9 shows the Nyquist plots of the diagonal open-loop processes of L(s) with the Gershgorin circles.The small circles in l11(s) and l22(s) indicate good diagonal dominance in these loops, not so in l33(s) with very large circles, which causes this loop to exhibit significant interaction.Figure 10 shows the closed-loop time responses achieved with the proposed MIMO PID control and the controllers of the two aforementioned authors.In the simulation, there is a unit step change in references 1, 2, and 3 at t = 1 min, t = 100 min, and t = 200 min, respectively.In addition, there is a 0.5 step change in all the process inputs at t = 300 min as a load disturbance.Table 2 lists the corresponding IAE and TV values.In the reference tracking performance in loops 1 and 2, the proposed control has a slightly slower rise time than the other loops, but similar settling time and no overshoot.In loop 3, the proposed control and Ram's controller are faster than Silva's control.In terms of disturbance rejection, the three controllers have a similar response.In terms of decoupling performance, the proposed control outperforms the other two, because owing to the two decoupling constraints imposed in the optimization, it has the best decoupling performance.Because of this, the proposed control obtains the lowest IAE values in loops 2 and 3, where the interactions are stronger.The IAE values in loops 2 and 3 of the proposed control are 30.6%and 4.6% lower than those of Silva, respectively, and 5.9% and 83% lower than those of Ram's control.However, the IAE value in loop 1 is 4.6% higher than that of Silva and 29.5% higher than that of Ram.The proposed control has the lowest TV values in loops 1 and 2, which are, respectively, 5% and 2.3% lower than those of Silva and 14.2% and 34.1% lower than those of Ram.However, it achieves the highest TV value in loop 3, which is 14.4% higher than that of Silva and 42% higher than that of Ram.

Example 3: 2 × 3 Shell Process
In this example, the system to be controlled is a non-square subset of the Shell process with two outputs and three inputs representing a heavy oil fractionator [36].This is provided in the 2 × 3 matrix GSH(s) in (30).The process, besides being a multivariable nonsquare system, shows significant interaction and contains significant time delays.The time constants are expressed in minutes.A frequency response array of 1000 elements in the range of [10 −6 , 1] rad/min is obtained from this matrix.The IAE values in loops 2 and 3 of the proposed control are 30.6%and 4.6% lower than those of Silva, respectively, and 5.9% and 83% lower than those of Ram's control.However, the IAE value in loop 1 is 4.6% higher than that of Silva and 29.5% higher than that of Ram.The proposed control has the lowest TV values in loops 1 and 2, which are, respectively, 5% and 2.3% lower than those of Silva and 14.2% and 34.1% lower than those of Ram.However, it achieves the highest TV value in loop 3, which is 14.4% higher than that of Silva and 42% higher than that of Ram.

Example 3: 2 × 3 Shell Process
In this example, the system to be controlled is a non-square subset of the Shell process with two outputs and three inputs representing a heavy oil fractionator [36].This is provided in the 2 × 3 matrix G SH (s) in (30).The process, besides being a multivariable non-square system, shows significant interaction and contains significant time delays.The time constants are expressed in minutes.A frequency response array of 1000 elements in the range of [10 −6 , 1] rad/min is obtained from this matrix.
To design the multivariable PID control, optimization II is used with bandwidth design specifications of 0.0075 and 0.012 rad/min for loops 1 and 2, respectively, in order to obtain cut-off frequencies similar to those of other authors.In addition, values of α = 85 • and β = 15 • are specified for both loops, and the two decoupling constraints (d) and (e) shown in (23) are used.
1.77e −84s 60s+1 5.88e −81s 50s+1 5.39e −54s 50s+1 5.72e −42s 60s+1 6.9e −45s 40s+1 (30) Taking the pseudoinverse of G SH (0) as the initial control K 0 , the proposed method obtains the 3 × 2 PID control shown in (31)    The Nyquist plots of the two resulting ELTFs with the r and q lines associated with the design constraints are shown at the top of Figure 12.Table 3 lists the other classical robustness margins and crossover frequency ωc achieved with the proposed design with those obtained by two other authors.Specifically, the table also shows the values obtained with the centralized PI control of Ghosh and Pan [36], and the multiloop PID controller of Jin et al. [32].The proposed control achieves similar phase margins, better gain margins, and better maximum sensitivity values.It also achieves the highest bandwidth in loop 1. Figure 12 shows the Nyquist plots of the diagonal elements of the resulting open-loop matrix L(s) with the corresponding Gershgorin circles.The circles with near-zero radii indicate that diagonal dominance is achieved in the system.The Nyquist plots of the two resulting ELTFs with the r and q lines associated with the design constraints are shown at the top of Figure 12.Table 3 lists the other classical robustness margins and crossover frequency ω c achieved with the proposed design with those obtained by two other authors.Specifically, the table also shows the values obtained with the centralized PI control of Ghosh and Pan [36], and the multiloop PID controller of Jin et al. [32].The proposed control achieves similar phase margins, better gain margins, and better maximum sensitivity values.It also achieves the highest bandwidth in loop 1. Figure 12 shows the Nyquist plots of the diagonal elements of the resulting open-loop matrix L(s) with the corresponding Gershgorin circles.The circles with near-zero radii indicate that diagonal dominance is achieved in the system.The closed-loop time response of the proposed 3 × 2 PID control compared with the other two above-mentioned designs is shown in Figure 13.The simulations were performed using a unit step change in the references at t = 0.1 min in the first reference and at t = 1000 min in the second reference.In addition, at t = 2000 min, there is a load disturbance consisting of a 0.1 step change at the three process inputs.The proposed control and Ghosh's centralized PI control have a similar response with similar rise times and settling times and excellent decoupling performance; the proposed control is slightly better in disturbance rejection.Jin's control shows more interaction, which is logical as it is a multiloop PID control.From the simulations, we calculated the IAE and TV values, which are shown in Table 3.In the TV column, there are three values for each design because there are three control signals.The proposed control achieves the lowest IAE values, which are 9.5% and 3.9% lower in loops 1 and 2 than those of Ghosh, respectively, and 13.3% and 43.8 lower than those of Jin.However, it has higher TV values.The TV values obtained with the proposed control in control signals 1 and 2, respectively, are 187% and 21% higher than those of Ghosh; the TV of the third signal is 39% lower than that of Ghosh.Compared with Jin's control, these TV values are 35%, 8.6%, and 243% higher in control signals 1, 2, and 3, respectively.Because the proposed methodology does not consider the input sensitivity function, the resulting controller can generate large control signals when the closed-loop bandwidth is chosen to be much larger than the open-loop bandwidth [40].The closed-loop time response of the proposed 3 × 2 PID control compared with the other two above-mentioned designs is shown in Figure 13.The simulations were performed using a unit step change in the references at t = 0.1 min in the first reference and at t = 1000 min in the second reference.In addition, at t = 2000 min, there is a load disturbance consisting of a 0.1 step change at the three process inputs.The proposed control and Ghosh's centralized PI control have a similar response with similar rise times and settling times and excellent decoupling performance; the proposed control is slightly better in disturbance rejection.Jin's control shows more interaction, which is logical as it is a multiloop PID control.From the simulations, we calculated the IAE and TV values, which are shown in Table 3.In the TV column, there are three values for each design because there are three control signals.The proposed control achieves the lowest IAE values, which are 9.5% and 3.9% lower in loops 1 and 2 than those of Ghosh, respectively, and 13.3% and 43.8 lower than those of Jin.However, it has higher TV values.The TV values obtained with the proposed control in control signals 1 and 2, respectively, are 187% and 21% higher than those of Ghosh; the TV of the third signal is 39% lower than that of Ghosh.Compared with Jin's control, these TV values are 35%, 8.6%, and 243% higher in control signals 1, 2, and 3, respectively.Because the proposed methodology does not consider the input sensitivity function, the resulting controller can generate large control signals when the closed-loop bandwidth is chosen to be much larger than the open-loop bandwidth [40].

Example 4: Experimental Quadruple Tank System
In this case, the two proposed methods (I and II) are applied to an experimental process: the quadruple tank plant [46] of the process control lab of the University of Cordoba.It is a 2 × 2 system where the controlled variables are the levels, h1 and h2 in cm, of the two lower tanks, and the control signals are the flow references, qL and qR in cm 3 /s, of the secondary loops that command the pumps.The plant is configured to have interaction problems and a multivariable right-half plane (RHP) of zero.The transfer matrix in (32) shows the identified model of the process around the operation point h = [17,18] cm and q = [135, 135] cm 3 /s.The time constants are expressed in seconds.The model in (32) From this matrix, a frequency response array of 1000 elements in the range of [10 −6 , 1] rad/s is obtained.To design the multivariable PID control with the proposed optimization I, the robustness specification provided with an Lm value of 0.8 and α value of 75° is used for both loops.In addition to the static decoupling, the decoupling constraint d is applied at a frequency of 0.002 rad/s in each loop.To apply the proposed optimization II, a bandwidth design specification of 0.004 rad/s is used for both loops.In addition, values of α = 75° and β = 35° are specified for both loops, and the two decoupling constraints (d and e) shown in (23) are used.Taking the inverse of GT(0) as the initial control K 0 , the proposed methods I and II obtain the 2 × 2 PID control shown in (33) and (34), respectively.In the case of optimization II, the linear margin values achieved are ℓ1 = 0.64 and ℓ2 = 0.63.For comparison, the previous iterative methodology by Garrido et al. in [39], which is extended in this study, is also applied to this example with a specification provided with

Example 4: Experimental Quadruple Tank System
In this case, the two proposed methods (I and II) are applied to an experimental process: the quadruple tank [46] of the process control lab of the University of Cordoba.It is a 2 × 2 system where the controlled variables are the levels, h 1 and h 2 in cm, of the two lower tanks, and the control signals are the flow references, q L and q R in cm 3 /s, of the secondary loops that command the pumps.The plant is configured to have interaction problems and a multivariable right-half plane (RHP) of zero.The transfer matrix in (32) shows the identified model of the process around the operation point h = [17,18] From this matrix, a frequency response array of 1000 elements in the range of [10 −6 , 1] rad/s is obtained.To design the multivariable PID control with the proposed optimization I, the robustness specification provided with an L m value of 0.8 and α value of 75 • is used for both loops.In addition to the static decoupling, the decoupling constraint d is applied at a frequency of 0.002 rad/s in each loop.To apply the proposed optimization II, a bandwidth design specification of 0.004 rad/s is used for both loops.In addition, values of α = 75 • and β = 35 • are specified for both loops, and the two decoupling constraints (d and e) shown in (23) are used.Taking the inverse of G T (0) as the initial control K 0 , the proposed methods I and II obtain the 2 × 2 PID control shown in (33) and (34), respectively.In the case of optimization II, the linear margin values achieved are ℓ 1 = 0.64 and ℓ 2 = 0.63.For comparison, the previous iterative methodology by Garrido et al. in [39], which is extended in this study, is also applied to this example with a specification provided with L m = 0.8 and α = 75 • in both loops.The resultant MIMO PID control K G (s) is shown in (35 Figure 14 shows the progress of the PID parameters during the iterative process in the three methodologies.The fact of performing the design of all the PID control parameters of K(s) in each iteration and the new decoupling constraints help the convergence of the proposed methods.The methods performed with optimization I (proposed I) and optimization II (proposed II) require eight and nine iterations, respectively, whereas Garrido's previous method converges in eleven iterations.Table 4 lists the classical robustness margins and crossover frequencies ω c achieved with the three designs.Proposed I and Garrido's design obtain almost the same values of these margins because both are very similar and aim to maximize the integral gains while fulfilling the same robustness margin specification.However, the proposed I design shows lower bandwidth frequencies than Garrido's; this is a consequence of achieving better diagonal dominance and better decoupling performance imposed by the new decoupling constraints in optimization I.The proposed II design is performed by specifying a higher bandwidth frequency than that obtained with optimization I.As a result, the proposed II control obtains the highest bandwidth frequencies but achieves lower robustness margins than the other two designs.The closed-loop time responses of the multivariable PID designs are shown in Figure 15.The experiments were performed using a 4 cm step change at t = 50 s in the first reference and at t = 3600 s in the second reference.The inverse response to the step change in the references is due to multivariable RHP zero.The proposed I design has the slowest rise times but the best decoupling performance.The proposed II design shows the fastest response in both loops.Garrido's control shows a very similar response in loop 2 as the proposed II design; however, in loop 1, it has a slower response than the proposed II design and shows more interaction than the other two designs.

Conclusions
This paper proposes a new iterative method for centralized PID control design that can be applied to stable processes.The methodology is based on the loop shaping of the Nyquist plot of the equivalent loop transfer function (ELTF) in centralized control, which captures the dynamics and interactions with the other loops.By working with arrays of The IAE and TV values calculated from the experimental data are listed in Table 4.The proposed II control presents the lowest IAE values, which are 27% and 9.7% lower than those of the proposed I design in loops 1 and 2, respectively, and 39.7% and 4% lower than those of Garrido's design.On the other hand, the proposed I design has the lowest TV values, which are approximately 32-34% lower than those of the proposed II control and those of the Garrido'22 design.

Conclusions
This paper proposes a new iterative method for centralized PID control design that can be applied to stable processes.The methodology is based on the loop shaping of the Nyquist plot of the equivalent loop transfer function (ELTF) in centralized control, which captures the dynamics and interactions with the other loops.By working with arrays of the frequency response, the method does not require approximations, even when the process has time delays or irrational functions.The method is an extension of our previous methodology in [39].Its main improvements are the following: • The complete design of the entire control matrix is formulated at each iteration as a single linear programming problem in which the frequency response of the ELTFs and the open-loop transfer functions are linearly parameterized with respect to the PID parameter vector.In this way, a better convergence is sought.In our previous work, each column of the controller matrix was separately designed; • Two constraints are proposed to improve the decoupling performance: static decou- pling at ω = 0 and decoupling at a frequency close to the crossover frequency of the ELTFs.These constraints not only manage to considerably reduce the interactions but also help to improve the convergence of the iterative procedure; • Two possible optimizations are proposed according to the design specifications: (I) maximize the integral gains to optimize the disturbance rejection response while complying with a minimum linear robustness margin in each loop and (II) maximize the linear robustness margins of each loop while ensuring minimum bandwidth in each loop.In our previous work, the only option was optimization I; • The new methodology has been extended not only for square processes but also for non-square systems.
The effectiveness of the proposed method is demonstrated in simulations on three processes of 2 × 2, 3 × 3, and 2 × 3 dimensions as well as a real quadruple tank system.Comparison with other authors' designs indicated that the proposed methodology achieves similar or superior performance.The main limitations of the methodology are as follows: (a) convergence is not completely assured because the design specifications may not be achievable; (b) the method is restricted to stable processes because the linear margin is not valid for unstable systems; and (c) the magnitude of the control signals is not considered in the optimization, so the designed control may produce high control signals when the bandwidth in the closed-loop is too high.

Figure 1 .
Figure 1.General scheme of the centralized control system.Figure 1.General scheme of the centralized control system.

Figure 1 .
Figure 1.General scheme of the centralized control system.Figure 1.General scheme of the centralized control system.

3 .
Recalculate the new open-loop matrix frequency response L m (jω) with the new PID MIMO control matrix.With these open-loop transfer functions, we calculate the new ELTFs and check if the design requirements are met in all the loops; 4.

Figure 4 .
Figure 4. Flow chart of the proposed iterative procedure.

Figure 4 .
Figure 4. Flow chart of the proposed iterative procedure.

Figure 5 .
Figure 5. Progression of the PID parameters and linear margins through the iterations in example 1.

Figure 5 .
Figure 5. Progression of the PID parameters and linear margins through the iterations in example 1.

Figure 6 .
Figure 6.Nyquist diagrams (in blue color) of the ELTFs in example 1 at the top and of the diagonal open-loop processes with their respective Gershgorin circles (in red color) at the bottom.

Figure 6 .
Figure 6.Nyquist diagrams (in blue color) of the ELTFs in example 1 at the top and of the diagonal open-loop processes with their respective Gershgorin circles (in red color) at the bottom.

Figure 7 .
Figure 7. Controlled outputs and control signals corresponding to the time responses in example 1.

Figure 7 .
Figure 7. Controlled outputs and control signals corresponding to the time responses in example 1.

Figure 8 .
Figure 8. Progression of the PID parameters through the iterations in example 2.

Figure 8 .
Figure 8. Progression of the PID through the iterations in example 2.

Figure 9 .
Figure 9. Nyquist diagrams (in blue color) of the ELTFs in example 2 at the top and of the diagonal open-loop processes with their respective Gershgorin circles (in red color) at the bottom.

Figure 9 .
Figure 9. Nyquist diagrams (in blue color) of the ELTFs in example 2 at the top and of the diagonal open-loop processes with their respective Gershgorin circles (in red color) at the bottom.

Figure 10 .
Figure 10.Controlled outputs and control signals for time responses in example 2.

Figure 10 .
Figure 10.Controlled outputs and control signals for time responses in example 2.
figure shows the progress of the PID parameters and margins Lm during the iterative process.The linear margin values achieved are ℓ1 = 0.666 and ℓ2 = 0.728.

Figure 11 .
Figure 11.Progression of the PID parameters through the iterations in example 3.

Figure 11 .
Figure 11.Progression of the PID parameters through the iterations in example 3.

Electronics 2024 , 27 Figure 12 .
Figure 12.Nyquist diagrams (in blue color) of the ELTFs in example 3 at the top and of the diagonal open-loop processes with their respective Gershgorin circles (in red color) at the bottom.

Figure 12 .
Figure 12.Nyquist diagrams (in blue color) of the ELTFs in example 3 at the top and of the diagonal open-loop processes with their respective Gershgorin circles (in red color) at the bottom.

Electronics 2024 , 27 Figure 13 .
Figure 13.Controlled outputs and control signals for the time responses in example 3.

Figure 13 .
Figure 13.Controlled outputs and control signals for the time responses in example 3.
cm and q = [135, 135] cm 3 /s.The time constants are expressed in seconds.The model in (32) has an RGA with diagonal elements equal to −0.21 and a multivariable RHP zero at s = 1/164.7.

27 Figure 14 .
Figure 14.Progression of the PID parameters through the iterations in example 4. The closed-loop time responses of the multivariable PID designs are shown in Figure 15.The experiments were performed using a 4 cm step change at t = 50 s in the first reference and at t = 3600 s in the second reference.The inverse response to the step change in the references is due to multivariable RHP zero.The proposed I design has the slowest

Figure 14 .
Figure 14.Progression of the PID parameters through the iterations in example 4.

Electronics 2024 , 27 Figure 15 .
Figure 15.Controlled outputs and control signals for the time responses in example 4.

Figure 15 .
Figure 15.Controlled outputs and control signals for the time responses in example 4.

Table 1 .
Robustness and performance indices in example 1.

Table 1 .
Robustness and performance indices in example 1.

Table 2 .
Robustness and performance indices in example 2.

Table 2 .
Robustness and performance indices in example 2.

Table 3 .
Robustness and performance indices in example 3.

s ω cp (rad/min) IAE TV
has an RGA with diagonal elements equal to −0.21 and a multivariable RHP zero at s = 1/164.7. ).

Table 4 .
Robustness and performance indices in example 4.