Iterative Method for Tuning Multiloop PID Controllers Based on Single Loop Robustness Speciﬁcations in the Frequency Domain

: Multiloop proportional-integral-derivative (PID) controllers are widely used for controlling multivariable processes due to their understandability, simplicity and other practical advantages. The main difﬁculty of the methodologies using this approach is the fact that the controllers of different loops interact each other. Thus, the knowledge of the controllers in the other loops is necessary for the evaluation of one loop. This work proposes an iterative design methodology of multiloop PID controllers for stable multivariable systems. The controllers in each step are tuned using single-input single-output (SISO) methods for the corresponding effective open loop process (EOP), which considers the interaction of the other loops closed with the controllers of the previous step. The methodology uses a frequency response matrix representation of the system to avoid process approximations in the case of elements with time delays or complicated EOPs. Consequently, different robustness margins on the frequency domain are proposed as speciﬁcations: phase margin, gain margin, phase and gain margin combination, sensitivity margin and linear margin. For each case, a PID tuning method is described and detailed for the iterative methodology. The proposals are exempliﬁed with two simulations systems where the obtained performance is similar or better than that achieved by other authors.


Introduction
Traditional proportional-integral-derivative (PID) controllers have been the preferred control algorithm in most industrial applications for decades because of their acceptable control effect, satisfactory robustness and simple structure [1]. There are thousands of works about PID controller design methods and PID tuning rules, most of them developed for single input single output (SISO) loops. However, most industrial processes are multivariable systems consisting of multiple inputs and multiple outputs with couplings between them. These process interactions can seriously deteriorate the control system performance when they are not considered in the design stage.
There are mainly two control approaches for dealing with interactions in multivariable systems: decentralized (or multiloop) control, or centralized control. The multiloop approach first decides the different loops pairing the inputs and outputs according to some measurement such as the relative gain array (RGA) [2] with its many variants; then, a single controller is tuned for each loop (Figure 1). In contrast, using a centralized control, a full matrix controller is designed. Although this last approach can achieve a greater interaction reduction and a better performance than multiloop control, it should be implemented when it is justified since it shows some disadvantages over the first approach: a full matrix controller is designed. Although this last approach can achieve a greater interaction reduction and a better performance than multiloop control, it should be implemented when it is justified since it shows some disadvantages over the first approach: more complicated design, higher sensitivity to modeling errors and uncertainties, more difficulty of implementation and so on [3]. Consequently, the multiloop control is preferred when the interaction level and characteristics of the multivariable process allow one to obtain a good performance using it. An intermediate approach consists of designing a decoupling compensator to reduce the interactions of the original system and obtain a diagonal dominant apparent process and then calculating a multiloop control for this [4]. Multiloop PID controllers are widely utilized due to their simple control structure, robustness performance, easiness to handle loop failure, easy understandability, fewer tuning parameters and other practical advantages. There are different types procedures to carry out the design. Detuning methods design the PID controllers for each loop using SISO techniques and ignoring the interactions from other loops and then, detune the parameters until achieving some limit such as the biggest-log-modulus (BLT) [5], which determines the tradeoff between performance and stability. They are very simple; however, they do not exploit the multivariable information of the process model to obtain a better system response.
Sequential loop closing methods close the loops one after the other, usually starting with the fastest one. They take into account the process interaction of the loop in the closing of the next loop, and so on [6,7]. The final design depends on the order of loop closing and iteration procedures are necessary. In the iterative design methods, the full process model information is considered to tune the controllers and meet some SISO specification in each loop. Starting from an initial guess of controller parameters set, each controller is retuned independently with the other loops being closed using the controllers calculated in the previous step. The iterative procedure is carried out until the controller parameters converge or the loop specifications are achieved.
Among alternative methods, other methodologies perform the design without iterations by means of simultaneous equation solving [8,9] or nonlinear optimization procedures with matrix inequality constraints [10,11], evolutionary or bioinspired algorithms [12,13]. Other authors use independent design procedures such as the application of internal model control to the multiloop PID design [14,15], the combination of independent design and Nyquist stability analysis [16], the dominant pole placement in [17], the direct method in [18,19] based on effective open loop processes, the decomposition of the process into a sum of basic modes of the multiscale control scheme [20], the single-iteration strategy in [21] with gain and phase margin specifications and so on. The main advantage of these methods is that they do not need a prior knowledge of the controllers in the other loops and therefore the SISO controller is tuned independently by using some bounds to guarantee stability and performance. Nevertheless, they can show some disadvantages Multiloop PID controllers are widely utilized due to their simple control structure, robustness performance, easiness to handle loop failure, easy understandability, fewer tuning parameters and other practical advantages. There are different types procedures to carry out the design. Detuning methods design the PID controllers for each loop using SISO techniques and ignoring the interactions from other loops and then, detune the parameters until achieving some limit such as the biggest-log-modulus (BLT) [5], which determines the tradeoff between performance and stability. They are very simple; however, they do not exploit the multivariable information of the process model to obtain a better system response.
Sequential loop closing methods close the loops one after the other, usually starting with the fastest one. They take into account the process interaction of the loop in the closing of the next loop, and so on [6,7]. The final design depends on the order of loop closing and iteration procedures are necessary. In the iterative design methods, the full process model information is considered to tune the controllers and meet some SISO specification in each loop. Starting from an initial guess of controller parameters set, each controller is retuned independently with the other loops being closed using the controllers calculated in the previous step. The iterative procedure is carried out until the controller parameters converge or the loop specifications are achieved.
Among alternative methods, other methodologies perform the design without iterations by means of simultaneous equation solving [8,9] or nonlinear optimization procedures with matrix inequality constraints [10,11], evolutionary or bioinspired algorithms [12,13]. Other authors use independent design procedures such as the application of internal model control to the multiloop PID design [14,15], the combination of independent design and Nyquist stability analysis [16], the dominant pole placement in [17], the direct method in [18,19] based on effective open loop processes, the decomposition of the process into a sum of basic modes of the multiscale control scheme [20], the single-iteration strategy in [21] with gain and phase margin specifications and so on. The main advantage of these methods is that they do not need a prior knowledge of the controllers in the other loops and therefore the SISO controller is tuned independently by using some bounds to guarantee stability and performance. Nevertheless, they can show some disadvantages that can limit their application. Poor resulting performance can be derived due to not using the information of the controllers in the other loops to calculate the real equivalent SISO processes. In addition, several of these methods may require model reductions to be applied, even more in the case of systems with multiple time delays where the equivalent transfer functions of the process are irrational. Furthermore, they can achieve satisfactory performance for 2 × 2 processes; however, their extensions to higher dimensional systems are usually difficult because of complicated calculations or important model approximations.
The common difficulty that all these multiloop methodologies attempt to overcome is the fact the controllers interact with each other. Therefore, the evaluation of the real performance of one loop requires knowing the controller dynamics in other loops. This work proposes an iterative methodology for tuning multiloop PID controllers for stable multivariable processes assuming a previous proper input-output pairing. The method is based on a structural decomposition of the system into separate SISO loops. The controller of each individual loop is tuned according to the corresponding effective open loop process (EOP) that includes the interaction transmission of the other loops being closed and consequently, needs the set of controllers calculated in the previous iteration step. Therefore, the method considers the full information of the multivariable process model, which is based on a frequency response matrix representation to avoid approximations or reductions of the EOP in case of systems with multiple time delays. Due to this representation, various SISO robustness specifications on the frequency domain are proposed for tuning each independent loop using different SISO methodologies. An initial version of this methodology was introduced only for 2 × 2 processes in [22], where phase margin, gain margin or a combination of them were used as specifications for the PID tuning methodology in [23], and where the proposed designs obtained less conservative responses than those achieved by other frequency response methods based on Gershgorin's bands [24]. In this work, further research has been performed extending the method for n × n systems and considering new robustness specifications and SISO PID tuning methods.
The paper is structured as follows. Section 2 provides a preliminary background for the proposed methodology such as the effective open loop processes, the different robustness specifications covered in this work and their corresponding SISO PID tuning methods. In Section 3, the proposed methodology is described. Section 4 illustrates the proposal with several simulation examples and different robustness margins. Finally, conclusions are summarized in Section 5.

Effective Open Loop Processes
From Figure 1, the process is given by the n × n transfer matrix G(s) in (1), where g ij (s) represents the transfer function between output y i and input u j , and the decentralized controller is given by the diagonal matrix K(s) in (2), where k i (s) is the controller for the loop between output y i and input u i . The diagonal structure of K(s) is obtained assuming a diagonal input-output pairing for introducing the notation of the proposed methodology. However, the interaction analysis can recommend other input-output pairings. To maintain a diagonal K(s) structure in these cases, the G(s) process matrix must simply be reordered accordingly.
Processes 2021, 9, 140 4 of 25 Supposing that loop i is still open while the other n − 1 loops are closed using controllers tuned by some method, it is desired to tune the controller k i (s). To analyze the apparent open loop process between output y i and input u i under these conditions, the structural decomposition scheme in Figure 2 is proposed. In this scheme, K 1 (s) represents the controller k i (s) and K 2 (s) is the set of other n − 1 controllers in their corresponding closed loop. From this scheme, the real open loop transfer function between output y i and input u i , and named g i (s), can be obtained as follows in (3), where G 12 (s) is a row vector with the n − 1 elements of row i different from g ii (s), G 21 (s) is a column vector with n − 1 elements of column i different from g ii (s), G 22 (s) is a n − 1 × n − 1 matrix with the other elements of G(s), and I is the n − 1 order identity matrix. Equation (3) can also be expressed in terms of the G(s) elements, as shown in (4).
Processes 2021, 9, x FOR PEER REVIEW 4 of 28 Supposing that loop i is still open while the other n − 1 loops are closed using controllers tuned by some method, it is desired to tune the controller ki(s). To analyze the apparent open loop process between output yi and input ui under these conditions, the structural decomposition scheme in Figure 2    These g i (s) transfer functions are the effective open loop processes (EOPs) and they capture the effective dynamics of each open loop process while the other loops are closed with their corresponding controllers. They are a key concept of the proposed methodology because they allow one to decompose the control system into equivalent SISO control loops to be designed. Then, the stability of a decentralized control system can be achieved if the characteristic equation of each equivalent single loop, as shown in (5), does not contain zeros in the right half plane [25].
Unlike the diagonal processes g ii (s), the EOPs consider the interaction effects with the other loops. However, the g i (s) transfer functions are not an intrinsic property of the G(s) since they depend on the controllers tuned in the other loops. The controller design can be performed using SISO loops methods; nevertheless, the modification of any controller will affect the other loops, where it will be necessary to update the controller tuning. Hence the proposed methodology, and others, applies iterations for tuning the controllers until their parameters converge or the design specifications are achieved. The EOP expressions for 2 × 2 processes are given in (6) where their dependency on the alternate loop controller can be appreciated.

SISO PID Tuning Methods and Robutstness Specifications
There are many SISO methods to calculate a PID controller k i (s) for the corresponding EOP g i (s). However, according to (4), the EOPs of a multivariable process can be complicated transfer functions and given the diverse casuistry of models, it is impossible to use a particular PID tuning formula such as those for first order plus time delay systems, second order plus time delay systems and so on. Furthermore, in case of systems with multiple time delays, the resulting EOPs are usually irrational transfer functions. In these cases, some methods suggest obtaining simplified EOP transfer functions by approximation. In contrast, the proposed methodology works using an array of the frequency response of the EOPs, and hence no reduction is performed, and any arbitrary order system can be represented.
For this reason, the proposed SISO PID tuning methods to be used in the iterative design procedure must be based on this kind of representation and frequency domain specifications. In this work, the proposed PID tuning methodologies use the Nyquist diagram as a design tool assuming that the EOPs and controllers do not contain poles in the right half plane and therefore, stability is achieved if there are no encirclements to the critical point (−1, 0). The following robustness margins are proposed as specifications: phase margin, gain margin, maximum sensitivity, and linear margin. They measure in different ways "how far from instability" the nominal loop is. They are described below and are represented in Figure 3.
The maximum sensitivity Ms is the inverse of the shortest distance of the Nyquist plot to the critical point −1; therefore, it is related to the radius of the circle centered in −1 and tangent to the Nyquist curve (point C), which must be consequently outside this circle. This margin ensures lower bounds on ϕm and Am at the intersection points of this circle with the unit circle and the negative real axis, respectively. The value of Ms can be calcalculated according to (8) and typical values are in the range from 1.4 to 2 [26]. In [27], an alternative linear margin is proposed for robustness. It is defined by the line r that intersects the negative real axis at the point (−1 + ℓ, 0) with an angle α (as depicted in Figure 3). The Nyquist curve must be placed below this line for stability, guaranteeing the lower bounds on ϕm, Am and Ms given in (9). By means of this linear margin, the PID control design can be formulated as a linear programming problem.
Processes 2021, 9, 140 6 of 25 The maximum sensitivity M s is the inverse of the shortest distance of the Nyquist plot to the critical point −1; therefore, it is related to the radius of the circle centered in −1 and tangent to the Nyquist curve (point C), which must be consequently outside this circle. This margin ensures lower bounds on φ m and A m at the intersection points of this circle with the unit circle and the negative real axis, respectively. The value of M s can be calcalculated according to (8) and typical values are in the range from 1.4 to 2 [26]. In [27], an alternative linear margin is proposed for robustness. It is defined by the line r that intersects the negative real axis at the point (−1 + , 0) with an angle α (as depicted in Figure 3). The Nyquist curve must be placed below this line for stability, guaranteeing the lower bounds on φ m , A m and M s given in (9). By means of this linear margin, the PID control design can be formulated as a linear programming problem.
In the following subsections, two SISO PID tuning methods are briefly described as candidates to be used in the proposed decentralized methodology. The first one uses either φ m , A m or M s as specifications; the second one is intended for the linear margin.

PID Tuning by Moving a Point of the Nyquist Diagram
Given a process g(s) and a controller k(s), this method designs k(s) by moving a point of the Nyquist diagram of g(s), in the absence of the controller, to a desired target point on the Nyquist plot of the open loop transfer function L(s) = g(s)k(s). In other words, it forces the Nyquist diagram of L(s) to pass through a specific point of the complex plane at a design frequency. After defining this point and the frequency design ω d , the magnitude and argument of g(jω d ) and L(jω d ) are known. The magnitude and argument of k(s) and its complex expression at this frequency can be derived as follows: In [23], this procedure is used for the design of non-interactive PID controllers with phase margin and gain margin as specifications for stable processes. The transfer function of a non-interactive PID controller is given by (11), where K P is the proportional gain, T I is the integral time constant and T D , the derivative time constant.
Decomposing its frequency response at frequency ω d into real and imaginary parts and comparing with the last equation in (10), the conditions in (12) are obtained. Once the target point of L(jω) is defined in the Nyquist diagram, |k(ω d )| and ϕ k (ω d ) are obtained according to (10). Thus, the PID parameters can be calculated from (12). To ensure a specified phase margin φ m at frequency ω d , the target point at the Nyquist plot of L(jω) will be the point A in Figure 3. Then, the magnitude and argument of the controller are as follows: Similarly, to guarantee a given gain margin A m , the Nyquist diagram must pass through point B of Figure 3. The resultant conditions for the controller are given in (14).
In this methodology [23], the ratio β = T D /T I must also be specified. Substituting T D = βT I into the last equation of (12), a quadratic equation in T I is achieved. After solving this equation, the expressions to calculate the PID controller parameters are summarized as follows: Considering that a PID controller can add a phase in the range [−90 • , +90 • ], the range of possible design frequencies is limited for the phase difference between the desired target point on L(jω) and the process g(jω) according to the argument condition in (10). Any frequency ω k producing that −90 • ≤ ϕ k (ω k ) ≤ +90 • is a candidate for the design frequency. Between these possible frequencies, the method in [23] proposes choosing the frequency that allows it to maximize the integral gain (K I = K P /T I ). This integral gain is related to the closed loop performance in terms of load disturbance rejection, which can be expressed according to the integrated error. Maximizing K I leads to minimizing the integrated error. However, a solution is not guaranteed when using simultaneously phase and gain margins as a combined specification. In this case, to reach both margins with the same PID parameters, the process frequency response must fulfill a set of conditions described in [23] and the design frequency cannot usually be selected.
In [28], a similar procedure is proposed for the tuning PID controllers based on the sensitivity margin M s as robustness specification. In this case, the desired target point at frequency ω d is defined by the point C in Figure 3, that is, the tangency point of the Nyquist plot with the circle of radius 1/M s and centered at −1. Given the magnitude of this radius (1/M s) and the depicted angle θ with the tangency point, the coordinates of point C are given by (16), where ϕ C is the angle of L(jω d ) at this point with respect to the negative real axis.
From (16), the magnitude and phase of L(jω d ) can be calculated as shown in (17). Additionally, then, conditions in (18) for the magnitude and phase of the controller can be obtained.
Processes 2021, 9, 140 8 of 25 Similarly to the previous case, the PID controller parameters can be calculated using the Equation (15) for a given ratio β. In this case, it is also necessary to verify the tangency of the Nyquist diagram of L(jω) at the contact point C with the circle centered at −1. The tangency condition function TCF in (18) must be below a given tolerance value to obtain a proper design.
Again, any frequency ω k can be chosen as the frequency design as long as the phase ϕ k (ω k ) that must be added by the PID controller at this frequency is within the range [−90 • , +90 • ]. The method in [28] calculates a PID controller in each frequency of the aforementioned range and evaluates the tangent condition in each design. Among all valid designs, the one with maximum integral gain K I is selected for improving load disturbance rejection.
The particular cases of proportional-integral (PI) control (β = 0) and proportionalderivative (PD) control can also be derived from (12). In these cases, the corresponding controller time constant is obtained according to (20). For the PI control the design range frequency is limited to frequencies where the PI phase contribution In this last case, there is no integral gain and maximizing the proportional gain K P is used as a rule to select the design.
PI control :

PID Tuning by Linear Programming
This PID tuning method was proposed in [27] for stable systems. The method shapes the Nyquist diagram of the open loop process L(jω) forcing it to be placed below the line r in Figure 3. This line is defined by the intersection point in the negative real axis at −1 + , and the line slope given by angle α. It is also based on an array representation of the system frequency response. The design problem can be formulated as a linear programming optimization using a PID controller with the parallel structure in (21), where K D is the derivative gain (K D = K P T D ). This structure can be parameterized as in (22) where ρ is the controller parameter vector and ψ(s) is a vector depending only on s.
Then, any point at frequency ω k of the Nyquist plot of L(jω) can be expressed as a linear function of the controller parameters as follows: Processes 2021, 9, 140 9 of 25 Every point of L(jω) must be below this line r. Thus, there is a constraint at each frequency ω k of the L(jω) frequency response array considering the line r in Equation (24). When it is desired to maximize the integral gain K I in order to optimize the closed loop performance for load disturbance rejection, the optimization problem in (25) can be defined for given values of and α.

Proposed Methodology
This section describes the proposed iterative algorithm for the multiloop PID design. It supposes a decentralized control structure using a diagonal input-output pairing of the system elements. Therefore, a previous resolution of the pairing problem and an appropriate rearrangement of process matrix G(s) must be performed. The proposed method is aimed for n × n square stable processes. Its flow diagram is depicted in Figure 4; it is basically composed by the initialization of EOPs and the iteration loop through a set of steps: SISO PID controller design for current EOPs, recalculation of EOPs, analysis of cost index for specifications, and checking of specification achievement. In these steps, the superscript m of variables denotes the iteration number, while the subscript i indicates the loop number from 1 to n.
Processes 2021, 9, x FOR PEER REVIEW 10 of 28 system elements. Therefore, a previous resolution of the pairing problem and an appropriate rearrangement of process matrix G(s) must be performed. The proposed method is aimed for n × n square stable processes. Its flow diagram is depicted in Figure 4; it is basically composed by the initialization of EOPs and the iteration loop through a set of steps: SISO PID controller design for current EOPs, recalculation of EOPs, analysis of cost index for specifications, and checking of specification achievement. In these steps, the superscript m of variables denotes the iteration number, while the subscript i indicates the loop number from 1 to n.  First, the iteration variable m is set to zero and the initial EOP g 0 i (s) of each loop i must be obtained before entering the iteration loop. As there are no previous controllers in the first iteration, expressions (3) or (4) cannot be used for the first EOP calculation. In this work, it was proposed to use the reduced equivalent open loop process (REOP) in (25) as the initial EOP g 0 i (s), where controller terms are excluded [29]. Since PID controllers have integral action, assuming a closed loop perfect control, the REOPs can be approximated from the structural decomposition in (3) as follows: Then, a PID controller k m i (s) must be tuned for each EOP according to the required specifications by means of some SISO method. Next, the EOP of each loop is updated using these new controllers and expression (4). As the EOPs can be irrational transfer functions, the proposed methodology works using a frequency response array for EOP representation instead of a concrete transfer function. Consequently, the SISO PID tuning methodologies must allow one to perform the design from a frequential description of the process.
Using the new recalculated EOPs g m+1 i (jω) and the current k m i (jω) controllers, the achievement of desired specifications in all the loops must be evaluated and quantified into a cost index J m . If this index is below a user-defined tolerance, the specifications have been achieved and the current design is accepted; otherwise the procedure is iterated, the iteration variable is incremented and new PID controllers are tuned using the updated EOPs. The iteration process continues until a design is accepted, the cost index does not change in three consecutive iterations, or a maximum number of iterations is reached. In the last two cases, the design with the smallest J index is selected. The main drawback of the proposed procedure consists of not guaranteeing the convergence. However, it works properly in most tested examples.
The PID tuning method and the cost index definition depend on the required specification. In the next subsections, these issues are described in detail for each proposed robustness margin. It is also necessary to note that other SISO methods can be used as long as they can be performed using a process representation based on a frequency response array.

Phase Margin Case
In this case, for each loop i of a n × n process, a desired phase margin φ mi * and a controller time constant ratio β i must be provided as procedure inputs together with the process transfer matrix G(s). In any iteration, there are n EOPs g m i (jω) calculated from the previous iteration m − 1. The PID tuning procedure of each loop is based on the method in [23], which has been adapted to work with a frequency response array representation. It may be divided in the following steps:

1.
Given the frequency response array of the EOP g m i (jω), obtain the design frequency range for which the tuning of PID controller parameters is possible according to the argument condition in (10). The phase provided by the controller must be in the range from −90 to 90 • .

2.
Determine the magnitude and phase contributions of the controller according to (13) for each possible design frequency ω k . 3.
Calculate the PID parameters according to (15) for each frequency ω k . As a result, as many PID control candidates as design frequencies are obtained.

4.
Evaluate stability of each design and choose that one with the greatest integral gain.
After obtaining PID controllers for each g m i (jω), the EOPs are recalculated with them and the achieved phase margin φ mi is evaluated for each new pair k m i (jω) g m+1 i (jω). Then, the cost index defined in (27) is obtained. If its value is less than a given tolerance, the multiloop PID design is accepted. A tolerance of 0.015 by loop is used for this index in this work.

Gain Margin Case
This case is very similar to the previous one. A desired gain margin A mi * is specified for each loop instead of a phase margin. The magnitude and phase contributions of the con-troller are calculated by means of expressions (14). The cost index to validate the multiloop PID design is given by (28), where A mi is the obtained gain margin for k m i (jω) g m+1 i (jω).

Combined Phase and Gain Margin Case
When only phase margin (or gain margin) is specified as a single requirement, the achieved design cannot ensure the loop stability. This problem can be solved using as specification a combination of desired gain margin A mi * and phase margin φ mi * for each loop. If the phase margin is positive and the gain margin is greater than one, the loop will be stable. The procedure is very similar to previous cases with the same steps. Now the cost index for the decentralized PID control design is given by (29). It is composed of two terms for each loop: one for the phase margin and another for the gain margin. A double tolerance of the previous cases is used for this index by loop.
The main difference of this case is found in the SISO PID tuning algorithm to obtain simultaneously both margins in the loop. This achievement is not always possible, even more when these requirements are aimed in all the loops. Therefore, the proposed method relaxes the phase margin specification in one direction by means of a path-following algorithm. It is an adaptation of the analytical method in [23] using a frequency response array representation. The algorithm receives as inputs the desired gain margin A * mi and phase margin φ * mi for the loop i, the frequency response array of the EOP and the time constant ratio β for the PID controller of this loop. It provides as output the PID control k C i that minimizes the cost index (28) in loop i using the design frequency ω k . The flow diagram of the algorithm is shown in Figure 5 and it represents in detail the block of general methodology in Figure 4 where the PID controller k m i (s) is tuned for g m i (jω). Iteration superscripts m and loop subscript i have been omitted for clarity in Figure 5. The single PID tuning procedure consists of the following steps:

1.
First, the current relaxed phase margin specification φ C m is initialized to the original one φ * m . This phase margin φ C m will be relaxed if it is necessary. The gain margin specification will remain the same.

2.
Given the EOP g(jω), obtain the design range of frequencies ω k , which allows the solution to the PID control tuning for the current phase margin φ C m according to the argument condition in (10). This step is similar to that of the tuning procedure in the phase margin case. Note that φ * m and φ C m are still the same.

3.
Calculate a PID controller k k for each frequency ω k using φ C m as the specification. It is also analogous to the tuning procedure described in the phase margin case. A set of PID controllers is obtained with the size of the design frequencies array.

4.
Evaluate the cost index J k in (29) for each design k k analyzing k k (jω) g(jω). Select the PID design with the minimum cost index J k . This defines the current lowest index J C and its corresponding PID controller k C . 5.
If the cost index J C is below a user-defined tolerance, the original specifications A * m and φ * m are assumed to be achieved, the current PID control design is accepted, and the algorithm finishes without relaxing the phase margin specification. Otherwise it is necessary to enter in the main part of the flow diagram, which decides how to relax the phase margin requirement. 6.
The first time the procedure enters here, no information is available to decide whether the current relaxed margin specification φ C m must be relaxed by increment or decrement. Therefore, the steps from 2 to 4 are repeated for both higher and lower values of φ C m . A variation of ±1 is considered adequate as a tradeoff between accuracy and speed. Then, both branches of the main part of the flow diagram are executed simultaneously; in practice one branch is performed before the other being the order irrelevant. The left branch develops the PID designs for φ C m − 1 and evaluates each cost index J k using the original specifications A * m and φ * m . The minimum cost index is defined as J − and the associated PID controller as k − . Similarly, the right branch performs the PID designs for φ C m + 1 obtaining the best PID controller k + with the corresponding minimum cost index J + . 7.
Determine the search path by means of the three cost indices J − , J C and J + . There are three possibilities: a.
If J C is the minimum cost index, the best solution is achieved using the current relaxed phase margin φ C m with the PID controller k C , and the procedure finishes. b.
If J + is the minimum cost index, the design is improved by incrementing the current relaxed phase margin, so φ C m is incremented, and the procedure enters the loop of the right branch. As φ C m has been updated, some cost indices and their corresponding PID controllers must be also updated accordingly as follows: J − = J C , k − = k C , J C = J + and k C = k + . Next, the new index J + and controller k + are calculated following again the right branch steps. c.
If J − is the minimum cost index, the design is improved by decrementing the relaxed phase margin, and the algorithm follows the loop of the left branch. Consequently, cost indices and their corresponding PID controllers are updated as follows: J + = J C , k + = k C , J C = J − and k C = k − . Then, the new index J − and controller k − are obtained after the branch steps.

8.
The process iterates into the branch determined in the previous step until no improvement is obtained in the cost index. Then, J C results in the minimum cost index, its associated PID controller k C is selected as the final one and the procedure ends.

Sensitiviy Margin Case
The sensitivity margin M s is related to the minimum distance of the Nyquist plot to the critical point −1. This indicator is more confident for the relative stability than gain and phase margins. Although the gain margin and phase margin of the loop L(jω) are adequate, they do not imply necessarily that the Nyquist diagram is far enough from the critical point for a good robustness condition. Nevertheless, the sensitivity margin always provides a minimal value for the gain and phase margin. In this section, the sensitivity margin M s is used as specification in the proposed multiloop PID tuning procedure of Figure 4. In this case the cost index is given by (30). A tolerance of 0.025 by loop was used for this index in this work.
Again, the main difference of this case is found in the SISO PID tuning algorithm to achieve the specification in each loop. The method proposed in [28] provides a PID controller for a given M s and angle θ at the tangency point. However, a suitable design is not always ensured due to the tangency condition fulfillment in (19). For this reason, a modification of this method was proposed in this work. It combines the proposal of [28] and the relaxation method proposed in the previous subsection for the combination of phase and gain margins. In this case, the relaxation was performed over the tangency point angle θ. Although the PID tuning procedure needs to receive this angle as input, it is not a requirement for the robustness specification and therefore, it can be modified. The angle is limited to the range [0 • , 90 • ] with typical values between 0 and 45 • . The proposed tuning method relaxes this angle when no proper PID design can be found. An initial value of 25 • or 30 • is recommended.
The proposed method receives as inputs the desired sensitivity margin M * si and angle θ i for the loop i, the frequency response array of the EOP and the time constant ratio β for the PID controller of this loop. As outputs, the method provides the PID controller k C i that minimizes the cost index (30) in loop i using the design frequency ω k . It also provides the relaxed angle θ C i to be used in the next iteration. Figure 6 shows the flow diagram of the proposed algorithm. It is almost the same one depicted in Figure 5 with some differences and it is described in a similar way below:

1.
First, the current relaxed angle θ C is initialized to the original one θ. This angle θ C will be relaxed if needed and will be provided as an algorithm output.

2.
Given the g(jω), obtain the design range of frequencies ω k , which allows one to find a solution for the PID control tuning given the sensitivity margin M * s and angle θ C according to the argument condition in (10). For this, it is necessary to calculate the magnitude and phase of L(jω) at the tangency point by means of (16) and (17).

3.
Calculate a PID controller k k for each frequency ω k using M * s and θ C as specifications. First, the argument and phase to be provided by the controller are computed from (18); then, using expressions in (15), a set of PID controllers is obtained as large as the size of the design frequency array.

4.
Evaluate the tangency condition function TCF k in (19) for each k k (jω) g(jω). Among all designs with TFC k below a user-defined tolerance tolTC, the one with maximum integral gain K I for improving load disturbance rejection is selected. This defines the PID controller k C . The sensitivity margin M s of k C (jω) g(jω) is evaluated and the corresponding cost index J C is calculated using (30).

5.
If the cost index J C is below a user-defined tolerance, the specification M * s is achieved with the original angle θ, the current PID control design is acceptable, and the algorithm ends without relaxing the angle θ C . Otherwise it is necessary to enter in the main part of the flow diagram, which decides how to relax this angle. 6.
In similar way to the algorithm proposed in Section 3.3, the procedure needs more information to whether the relaxed angle θ C must be incremented or decremented. Therefore, steps from 2 to 4 are repeated for θ C + 1 and θ C − 1. The left branch performs the PID designs for θ C − 1 obtaining the PID controller k − and its corresponding cost index J − . The right branch develops the PID designs for θ C + 1 providing the PID controller k + with the corresponding cost index J + . 7.
Determine the search path by means of the three cost indices J − , J C and J + . There are three possibilities: a.
If J C is the minimum cost index, the best solution is achieved using the current relaxed angle θ C with the PID controller k C , and the procedure finishes. b.
If J + is the minimum cost index, the tuning is improved incrementing the relaxed angle, so θ C is incremented, and the process enters the loop of the right branch. Since θ C has been updated, some cost indices and their corresponding PID controllers must be also updated accordingly as follows: J − = J C , k − = k C , J C = J + and k C = k + . Next, the new index J + and controller k + are obtained following again the right branch steps. c.
If J − is the minimum cost index, a decrease in θ C improves the PID design, so it is decremented, and the algorithm follows the loop of the left branch. This θ C modification entails the update of the following indices and associated controllers: J + = J C , k + = k C , J C = J − and k C = k − . Then, the new index J − and controller k − are obtained after the left branch steps.

8.
The process iterates into the branch determined in the previous step until no improvement is achieved in the cost index. Then, J C results in the minimum cost index and the procedure ends. The algorithm provides as outputs the associated PID controller k C and the relaxed angle θ C to be used as input in the following iteration step of the general multiloop design procedure.
Processes 2021, 9, x FOR PEER REVIEW 16 of 28 Figure 6. Flow diagram of the proposed SISO PID tuning algorithm for a specified sensitivity margin.

Linear Margin Case
In this case, the monovariable PID tuning design is performed through a linear programming optimization as proposed in [27] for stable systems and commented in Section 2.2.2. The method uses a linear margin defined by the line r intersecting the negative real axis at the point (−1 + ℓ, 0) with an angle α and forces the Nyquist plot L(jω) to be located below this line. This ensures lower bounds on ϕm, Am and Ms. The application of this method into the proposed algorithm in Figure 4 for the multiloop PID design is straightforward. For each loop i of the n × n process G(s), the corresponding linear margin must be defined by ℓi and αi as inputs. Once the algorithm enters the iterative procedure of the

Linear Margin Case
In this case, the monovariable PID tuning design is performed through a linear programming optimization as proposed in [27] for stable systems and commented in Section 2.2.2. The method uses a linear margin defined by the line r intersecting the negative real axis at the point (−1 + , 0) with an angle α and forces the Nyquist plot L(jω) to be located below this line. This ensures lower bounds on φ m , A m and M s . The application of this method into the proposed algorithm in Figure 4 for the multiloop PID design is straightforward. For each loop i of the n × n process G(s), the corresponding linear margin must be defined by i and α i as inputs. Once the algorithm enters the iterative procedure of the flow diagram in Figure 4, the n EOPs g m i (jω) are already calculated from the previous iteration m − 1. Then, a PID controller is tuned for each EOP by means of the linear programming problem in (25). Next, the new EPO g m+1 i (jω) of each loop is updated using these new sets of PID controllers.
For each new pair k m i (jω) g m+1 i (jω) the linear margin must be evaluated to quantify a cost index to decide whether the iteration process continues or not. This is the main point that differs from previous cases. There are two possible situations: (a) all points of L(jω) are below line r or (b) some points are above it. If all L(jω) points are below the line r as represented in Figure 7a, it is preferred that the Nyquist plot to be as close as possible to the line r to avoid unnecessary conservative designs. This can be measured by the shortest distance d min of the Nyquist diagram to the line r, which is calculated selecting the minimum value of the distances of each point L(jω k ) to the line r. Considering the line r equation in (24), the distance of the point L(jω k ) to this line is given by (31). When there are some points of L(jω) above the line, as shown in Figure 7b, the linear margin is not fulfilled, and again it is desired that these points to be as close as possible to the line r. This can be measured by the greatest distance dmax of these Nyquist points to the line r, which is the maximum distance. Therefore, this work proposed to define the fitting distance D of L(jω) to line r as follows: Using this fitting distance to the line r, the cost index defined in (33) was proposed as the sum of the fitting distances of all the loops. If its value is below a user-defined tolerance (0.002 by loop), the multiloop PID tuning is approved and the algorithm ends; otherwise, a new iteration is performed. In [30], an iterative procedure for decentralized controllers is also proposed using the linear margin in [27]; however, it performs the iterations in other way and does not use a cost index for evaluating the design or the algorithm convergence. When there are some points of L(jω) above the line, as shown in Figure 7b, the linear margin is not fulfilled, and again it is desired that these points to be as close as possible to the line r. This can be measured by the greatest distance d max of these Nyquist points to the line r, which is the maximum distance. Therefore, this work proposed to define the fitting distance D of L(jω) to line r as follows: Using this fitting distance to the line r, the cost index defined in (33) was proposed as the sum of the fitting distances of all the loops. If its value is below a user-defined tolerance (0.002 by loop), the multiloop PID tuning is approved and the algorithm ends; otherwise, a new iteration is performed. In [30], an iterative procedure for decentralized controllers is also proposed using the linear margin in [27]; however, it performs the iterations in other way and does not use a cost index for evaluating the design or the algorithm convergence.

Illustrative Examples
In this section, the effectiveness of the proposed methodology for tuning multiloop PID controllers is demonstrated using two simulation processes. The multivariable systems are 2 × 2 and 3 × 3 stable processes. The methodology is applied using different SISO PID tuning methods.

Wood and Berry Distillation Column
The transfer function matrix of this 2 × 2 system with time delays is given by (34). It represents a binary distillation column process with top and bottom compositions as outputs to be controlled and reflux and steam flows as manipulated inputs [31]. This is a classical testbench in the multivariable control. The time constants and delays are expressed in minutes. The process RGA in the stationary state is shown in (35) and the diagonal positive elements recommend a diagonal input-output paring. A frequency response array of 1000 elements is calculated within the frequency range [10 −5 , 10] rad/s.
Five decentralized controllers are designed to verify the proposed methodology with the five SISO PID tuning methods described in Section 3. The required specifications in each case are the following ones: In each case, the PI controllers are firstly attempted to achieve the desired specifications. If they are not reached, it is checked whether the decentralized design is improved using PID controllers. A time constant ratio β = 0.1 is set for the PID control. Table 1 collects the PID controllers obtained in each case with the procedures proposed in this work. The provided PID parameters are those of the parallel PID structure in (21). The table also shows the obtained specifications of phase margin, gain margin, sensitivity margin and the phase margin crossover frequency, which is usually close to the loop bandwidth. PI controllers are enough for the phase margin case and sensitivity margin case. In all cases, the requirements are practically achieved in both loops. In the third design, the obtained phase margin for the second loop is 46 • instead of 45 • , but it is very close. Figure 8 displays the PID parameter values through the iteration procedure in each case. Correspondingly, Figure 9 shows how the required specifications were achieved and the cost index decreased and converged reaching a value below the user-defined tolerance. The designs were obtained between 5 and 10 iterations, being most of them very close to the final values after four iterations. The case of phase and gain margin combination shows the slowest convergence. The cases of phase margin and sensitivity margin were the fastest ones, which was logical since these cases used PI controllers and had fewer tuning parameters.   The resultant Nyquist diagrams of the proposed designs are illustrated in Figure 10, where each column corresponds to the two loops of each case. The closed loop system responses of the proposed decentralized controllers are shown in Figure 11. Unit step changes were applied at t = 1 min in the first reference and at t = 80 min, in the second one. There were also a 0.25 step change in both inputs as load disturbance at t = 160 min. In the performed simulations, the PID controllers were implemented using a PI-D structure, where the control signal was calculated according to (36) and the derivative action is filtered by a first order term with N = 20 [32]. The effect of this filter in the frequency response was neglected into the frequency range of interest.   The resultant Nyquist diagrams of the proposed designs are illustrated in Figure 10, where each column corresponds to the two loops of each case.  Figure 11. Unit step changes were applied at t = 1 min in the first reference and at t = 80 min, in the second one. There were also a 0.25 step change in both inputs as load disturbance at t = 160 min. In the performed simulations, the PID controllers were implemented using a PI-D structure, where the control signal was calculated according to (36) and the derivative action is filtered by a first order term with N = 20 [32]. The effect of this filter in the frequency response was neglected into the frequency range of interest. The integrated absolute errors (IAEs) of each simulation case and the total variation (TV) of the control signals were calculated. They are shown in the last two columns of Table 1. The worst two cases were the designs obtained only for the phase margin or only for the gain margin. They show a more oscillatory response and higher IAE and TV values. The multiloop PI controller achieved with only the phase margin as a requirement had a low value of gain margin in the second loop. In a similar way, the design performed only with gain margin as specification achieved a low phase margin in the first loop. This resulted in a poorer performance of the control system. The integrated absolute errors (IAEs) of each simulation case and the total variation (TV) of the control signals were calculated. They are shown in the last two columns of Table 1. The worst two cases were the designs obtained only for the phase margin or only for the gain margin. They show a more oscillatory response and higher IAE and TV values. The multiloop PI controller achieved with only the phase margin as a requirement had a low value of gain margin in the second loop. In a similar way, the design performed only This situation was improved using the combined specification of gain margin and phase margin, the sensitivity margin, or the linear margin. These three cases ensured lower bound values of the gain margin and phase margin. The proposed controllers of the last two cases were quite similar and they achieved the best performance with lower values of IAE and TV. The proposal obtained for the combination of φ m and A m had a good response for reference step changes; however, the bandwidth of its second loop was a bit low, which resulted in a slow load disturbance rejection and a high IAE value in the corresponding loop.

Ogunnaike and Ray Distillation Column
The transfer matrix of this 3 × 3 distillation column process is given by (37) in [33], where delays and time constants are expressed in minutes. According to its stationary RGA, which is shown in (38), a diagonal input-output paring is recommended since the diagonal elements were positive and close to the unity. From this model, a frequency response array of 1000 elements is obtained within the frequency range from 10 −5 to 10 rad/s. The proposed methodology was applied to this process using PI controllers and using two different specifications: firstly, with the sensitivity margin (proposed 1), and secondly, with the linear margin (proposed 2). The values specified for these requirements in each loop and the corresponding achieved robustness margins were collected in Table 2. Figure 12 shows the evolution of the PI parameters, the required specifications and the cost index value for each iteration. Proposal 1 converged after four iterations and the second one, after six. The resultant Nyquist diagrams of each loop are illustrated in Figure 13, where the first row corresponds to proposal 1 and the second one, to proposal 2.    In comparison with example 1, which has been used to illustrate the majority of aspects and cases of the proposed methodology, the proposals for example 2 were limited only to these two cases because the obtained multiloop PI controllers were compared with In comparison with example 1, which has been used to illustrate the majority of aspects and cases of the proposed methodology, the proposals for example 2 were limited only to these two cases because the obtained multiloop PI controllers were compared with two decentralized controllers proposed by other authors: the multiloop PI control of Vu in [19] and the multiloop PID control of Huang in [18]. Figures 14 and 15 show a comparison of the closed loop system responses where unit step changes were applied at t = 1 min in the first set-point, at t = 100 min, in the second one and at t = 200 min, in the third one. A unit step change was also applied in all inputs at t = 400 min as load disturbance. Figure 14 details the time response of the controlled variables while Figure 15 shows the control signals. Table 2 indicates the IAE values and TV values obtained by the four designs. It also includes the robustness margins obtained with the controllers proposed by the other authors.  The two proposed multiloop PI controllers had a better performance than that obtained by the decentralized PI control of Vu, which shows the highest TV values. The multiloop PID control of Huang achieved the lowest IAE values in loop 1 and 2 at the expense of a worse response of the loop 3 with a higher IAE value and a great interaction from loop 1. Processes 2021, 9, x FOR PEER REVIEW 26 of 28

Conclusions
An iterative multiloop PID tuning methodology for multivariable square and stable processes with multiple time delays was developed in this work. The iterative procedure consists of the following steps: a decomposition of the systems into separate SISO loops by means of the concept of effective open loop processes; the tuning of each PID controller for the corresponding EOP using a robustness margin in the frequency domain; and recalculation of EOPs with these new controllers and evaluation of the achieved specifications by a cost index. Due to the use of EOPs, the interaction from other loops and the full information of the multivariable process model are considered for tuning each PID controller. Since the EOP transfer functions can result in very complicated or irrational expressions, the proposed method is based on a frequency response matrix representation to avoid approximations or reductions of the EOPs. For this reason, the proposed SISO PID tuning methods to be used in the iterative design procedure are based on this kind of representation and robustness specifications in the frequency domain. The following cases were described and analyzed: • Phase margin or gain margin cases: the proposed SISO PID tuning method is based on moving a point of the EOP Nyquist diagram to a desired one of the open loop L(jω) fulfilling the specified margin and providing the maximum possible integral gain.

•
Combined phase margin and gain margin case: a new PID tuning method is proposed based on the previous cases. Since achieving both margins simultaneously is

Conclusions
An iterative multiloop PID tuning methodology for multivariable square and stable processes with multiple time delays was developed in this work. The iterative procedure consists of the following steps: a decomposition of the systems into separate SISO loops by means of the concept of effective open loop processes; the tuning of each PID controller for the corresponding EOP using a robustness margin in the frequency domain; and recalculation of EOPs with these new controllers and evaluation of the achieved specifications by a cost index. Due to the use of EOPs, the interaction from other loops and the full information of the multivariable process model are considered for tuning each PID controller. Since the EOP transfer functions can result in very complicated or irrational expressions, the proposed method is based on a frequency response matrix representation to avoid approximations or reductions of the EOPs. For this reason, the proposed SISO PID tuning methods to be used in the iterative design procedure are based on this kind of representation and robustness specifications in the frequency domain. The following cases were described and analyzed: by means of an iterative path-following algorithm until no improvement is achieved in a defined cost index that combines both margins. • Sensitivity margin case: a new PID tuning method is also proposed for this case. It is based again on moving a point of the Nyquist diagram to fulfill the specification. The algorithm receives as inputs the sensitivity margin M s and angle θ at the tangency point, which can be relaxed in similar way to the previous case. In the method this angle is relaxed if needed to achieve a suitable design fulfilling properly the tangency condition and maximizing the integral gain. • Linear margin case: the proposed PID tuning method performs a loop shaping of L(jω) approaching it as a linear optimization problem. It ensures that the Nyquist diagram of L(jω) is below a line r and maximizes the integral gain. In this work a new cost index based on the distance of the Nyquist plot to the line r is proposed to determine when the algorithm ends.
The different cases were illustrated through two simulations examples. The results were compared with other works and the proposed methodology obtained similar or better performance. The best results were achieved when the sensitivity margin or linear margin were used as specifications. In addition, they usually require less iterations for the design. The main disadvantage of the proposed iterative procedure is that convergence is not always ensured. In this work the iteration process ended when the cost index associated to the specifications was below a user-defined tolerance, the cost index did not change in three consecutive iterations, or a maximum number of iterations was reached. However, the proposed methods work properly in most tested processes achieving the required specifications or others close to them.