Iterative Design of Centralized PID Controllers Based on Equivalent Loop Transfer Functions and Linear Programming

This paper deals with the design problem of multivariable Proportional-Integral-Derivative (PID) controllers for square and stable multivariable processes when a linear margin at the Nyquist plot is considered as robustness specification for each closed loop. A tuning method is developed based on the new concept of equivalent loop transfer function, which is proposed for centralized control and allows an independent design for each loop considering the interactions with the other loops through an iterative procedure. For the k-th loop, the PID parameters of the k-th column of the control matrix are calculated in each iteration by a linear programming optimization that maximizes the integral gains while fulfilling the robustness specification and achieving static decoupling. The method uses a frequency response array as representation of the process, which allows its applicability to systems with multiple time delays without requiring model reductions or approximations. The effectiveness of the method is illustrated by means of two simulation examples with dimensions $2 \times 2$ and $3 \times 3$ . Comparisons with other centralized control methodologies show that the proposed approach achieves similar or greater performance and a remarkable better disturbance rejection response.


I. INTRODUCTION
Despite most industrial processes are multiple inputs and multiple outputs (MIMO) systems, most of the thousands of works about design methods of Proportional-Integral-Derivative (PID) control, which has been the preferred one in industrial process control for decades [1], have been developed only for single input single output (SISO) loops [2]. The MIMO control design becomes more complicated since it requires additional considerations to the SISO case, such as the process interactions between inputs and outputs or the presence of multiple time delays, to avoid a serious deterioration of the control system performance. Two control approaches are mostly used to cope with multivariable systems: decentralized control and centralized control.
The associate editor coordinating the review of this manuscript and approving it for publication was Fu-Kwun Wang . In the decentralized control scheme, the process inputs and outputs pairing problem is firstly solved and then, a SISO controller is designed for each loop. Consequently, the resultant decentralized control, also called multiloop control, is a diagonal matrix controller. This approach is satisfactory for weakly coupled processes. By contrast, for non-diagonally dominant systems with great couplings decentralized controllers have a worse performance as they cannot reduce enough the interaction level and important transitory responses can arise in other loops when some reference changes. To face such processes, a centralized control scheme, which uses a full matrix controller, is recommended. Centralized control approach can be developed in two ways: combining a decoupling compensator and a diagonal controller or directly designing a full matrix control. In the first one, the decoupler is calculated to reduce the process interactions and obtain a diagonally dominant apparent process; then, the elements of the diagonal controller are tuned with SISO methods [3]- [5]. Some authors combine both blocks to form a full matrix control [6]- [8]. In a pure centralized control scheme, the full matrix controller works as the only block to reduce interactions and control the different process outputs.
Literature provides methodologies for both decentralized and centralized control schemes. Some decentralized approaches are detuning methods [9], sequential loop closing method [10], optimization procedures [11]- [13], direct design [14], iterative methods [15], or intelligent control techniques [16]. These methods usually use PID elements and obtain a multiloop PID controller. For centralized control, some methodologies obtain non-PID elements in the controller matrix and are based on the paradigm of decoupling control [17], recursive least square optimization [18], H 2 optimal performance requirements [19], or H ∞ control procedures [20]. Some authors propose to approximate the resulting controllers to PID structures [21], unless they cannot fully meet the imposed requirements [22]. Other methods directly obtain a full matrix of PID controllers, that is a centralized PID controller, and are based on convex optimization [23], internal model control [24], or effective open loop transfer functions [25]- [28]. Although multivariable control design can also be developed from state space methodologies [29], they usually entail delay-free systems. Conversely, many industrial processes have time delays [30] and therefore, most of these centralized methods are based on a transfer matrix approach.
One problem of most of these multivariable methodologies is the increasing complexity of the elements for high dimensional systems, where irrational or complicated transfer functions can arise. They usually need model approximations or controller reductions that affect the robustness of the control system or can lead to conservative control designs. To avoid these approximations, some authors use the frequency response matrix representation of the system and robustness specifications on the frequency domain such as gain margin, phase margin, or modulus margin. These margins provide good robustness performance with respect to the uncertainties in both disturbance and process model, and are well accepted by process control engineers [31]. In [32], Karimi proposes a linear margin that ensures lower bounds for the previous classical margins, and develops an associated design methodology of SISO PID controllers based on linear programming. In [15], [33], this method is extended to decentralized PID controllers for stable MIMO systems by means of an iterative procedure.
The literature review shows that research on MIMO PID control design is still an important issue from an industrial application point of view. New design methods are desired to be simple, general, and robust. This work presents the new concept of equivalent loop transfer function (ELTF) of centralized control systems. The ELTF accounts for loop interactions and allows to design separately each controller matrix column in an iterative procedure when an initial guess of the controller matrix is provided. Based on the ELTF concept, a new iterative methodology for designing centralized PID controllers for stable and square MIMO systems is proposed and formulated for each loop as a linear programming optimization that maximizes integral gains of PID controllers subject to constraints on robustness linear margin and static decoupling. The proposed method has some advantages in comparison to other related methodologies: • Since a frequency response array is used as representation of the process, the method is applicable to general square and stable systems with non-minimum phase zeros, time delays or any order irrational transfer functions without requiring any approximation or model reduction.
• Lower bounds on robustness margins in the frequency domain can be specified for each ELTF.
• With the static decoupling constraint, the process interactions are reduced at low frequencies and the method can be applied to non-diagonally dominant systems.
• Compared to other centralized controllers, better disturbance rejection responses are obtained because of the integral gain maximization. The outline of the paper is as follows: section II describes the concept of ELTF and its linear parametrization, the problem formulation, and the proposed iterative procedure. Although it is mainly expounded for 2 × 2 processes, expressions for higher dimensional systems are also provided. Simulation results in comparison with other methods are illustrated in Section III to verify the effectiveness of the proposed method. Conclusions are summarized in Section IV.  Fig. 1 shows the unity feedback control scheme for a centralized control system, where the n × n transfer matrix G(s) represents the MIMO process, the n × n transfer matrix K(s) is the full matrix controller, the process element g ij (s) is the transfer function between output y i and input u j , the transfer function k ij (s) is the control element between output y j and control signal u i , and the signal r i is the reference for process output y i . The difference between r i and y i is the error e i . The loop transfer matrix L(s) is defined as the transfer matrix around the loops as seen from the outputs when breaking all the loops. When unity feedback is assumed, L(s) = G(s)·K(s). Then, for the 2 × 2 case, the L(s) elements are calculated according to (1). In Fig. 2, the proposed equivalent loop transfer function l 1 (s) for the first loop is shown. This ELTF represents the transfer function between output y 1 and error signal e 1 when the first loop is open and the other one is closed. Unlike loop transfer functions l ii (s) in (1), the ELTF captures the effective dynamics of the corresponding loop transfer function and interactions with other loops. The concept of ETLF can be considered an extension for centralized control of the equivalent open-loop process for decentralized control in [15], [33]. For the 2 × 2 case, the expressions of the ELTFs l 1 (s) and l 2 (s) are given by (2).

L(s)
From (1), the elements l ij (s) of the j-th column of L(s) depend on the controllers k ij (s) of the j-th column of K(s). This dependency by columns allows to represent the ELTF l j (s) regarding the elements k ij (s) of the j-th column of K(s) when the other columns of K(s) are assumed to be known (and consequently, the other non j-th columns of L(s) are also known). For instance, by substituting the elements l 11 (s) and l 21 (s) given in (1) into the expression of l 1 (s) in (2) we obtain the new equation for l 1 (s) in (3). Therefore, the ELTF l 1 (s) depends only on k 11 (s) and k 21 (s) if l 12 (s) and l 22 (s) are assumed to be known. Then, k 11 (s) and k 21 (s), which are the controller elements of the first column of K(s), must be designed to achieve the desired loop specifications of the first loop as well as the stability of the control system. The stability is reached if right half plane zeros are avoided in each closed loop characteristic equation: 1 + l i (s) = 0.
In the general case of n × n processes, the elements l ij (s) of L(s) are given by (4). The ELTF l j (s) is calculated according to (5), which shows that the ETLF l j (s) can be expressed again involving only the controller elements of the j-th column of K(s) when the other columns of K(s), and therefore also of L(s), are known.
The design of the controller elements of a same column must be performed simultaneously to achieve the desired specifications of the corresponding loop. It can be developed independently of the other loops and columns of K(s). However, the modification of any column of K(s) will affect the other loops, and consequently, the tuning of the corresponding controllers must be updated. Therefore, the proposed methodology performs an iterative tuning procedure until the design requirements are achieved in all loops, as explained later.
To obtain a centralized PID controller, the elements of the control matrix K(s) are proposed with the parallel PID structure given in (6), where K Pij , K Iij and K Dij are the proportional gain, integral gain, and derivative gain, respectively, of k ij (s). With the parallel structure, every control element k ij (s) can be parameterized as shown in (7), where ρ ij denotes the vector of control parameters of k ij (s) and ψ(s) is a vector that depends only on s. Then, the product of the process element g km (s) and the controller element k ij (s) is parameterized as follows in (8).
For the case 2 × 2, the PID parameter vectors of each column of K(s) are defined by ρ 1 and ρ 2 in (9). By substituting (8) into (1) and using (9) we express L(s) according to (10), where ϕ ij (s) is a vector of six elements that depends only on s and is used to shorten notation. Note that, as shown in (10), the elements of the j-th column of L(s) are parameterized on ρ j .
Combining (10) with (2) gives a linear parameterization of the ELTFs l 1 (s) and l 2 (s) as shown in (11). For instance, in l 1 (s) the elements l 11 (s) and l 21 (s) in (2) are interchanged by those in (10) whereas l 12 (s) and l 22 (s) are assumed to be previously known. In the case of l 2 (s), the elements of the second column of L(s) are replaced and those of the first column are considered known. Note that the ELTF l j (s) depends only on the parameter vector ρ j , which is associated to the j-th column of K(s), while χ j (s) is a vector of six elements that depends only on information from the process G(s) and a previous matrix L(s), which implies a previous control matrix K(s).
For higher dimensional systems, a similar procedure can be performed by expressions (4) and (5) to obtain the elements of L(s) and the ELTFs as linear parameterizations in the vectors of PID parameters of each column of K(s). For instance, in the case of 3 × 3 processes, the elements l 11 (s) and the ELTF of loop 1 are parameterized as shown in (12) and (13), respectively, where the vector ρ 1 is defined as ( ρ T 11 ρ T 21 ρ T 31 ) T and has nine parameters, ϕ i1 (s) and χ 1 (s) are vectors of nine elements, and the loop transfer functions l 12 (s), l 13 (s), l 22 (s) and l 33 (s) are assumed to be known.
According to (2) and (5), the ELTF l i (s) of a MIMO process can result in a very complicated transfer function, even irrational as happens when there are time delays. For some design methodologies, this fact may require model reductions or complicated tuning procedures. In view of this fact, the proposed method uses a frequency response array of G(s) in the frequency range of interest, and thus, it can be applied to any arbitrary order system or time delay processes without requiring approximations. Therefore, every point at frequency ω k of the Nyquist diagram of the ETLF l j (jω k ) is stated as a linear function of the vector ρ j and is decomposed into real and imaginary parts, as shown in (14) for l 1 (jω k ) in the 2 × 2 case. A similar procedure is performed for n × n systems.

B. OPTIMIZATION BY LINEAR PROGRAMMING
Since the proposed methodology is based on a frequency response representation, frequency domain specifications on the Nyquist diagram are preferable as design requirements, and specifically, the robust linear margin described in [32] is selected in this work. This margin is determined by the line r that intersects the negative real axis at the point (−1 + , 0) with an angle α, as shown in Fig. 3. For stable processes, closed loop stability is guaranteed if the critical point (−1, 0) is not encircled; therefore, the Nyquist diagram must be positioned below this line providing stability and lower bounds on the phase margin φ m , gain margin A m , and modulus margin M m , which are represented by points A, B, and C, respectively. The lower bounds provided by the  [32]. If lower bounds for two conventional margins are specified, the values of α and that ensure such bounds can be calculated from (15) or from Fig. 4. As the values of α or increases, the system becomes better damped with less overshoot, less control effort, and better robustness; however, generally the bandwidth decreases and the settling time increases.
When this linear margin is used as design specification for the ELTF l j (jω), the PID controller parameters of the jth column of K(s), that is, the vector ρ j , can be calculated by a linear programming optimization if initial guess values for the other columns of K(s) are defined, and therefore, the corresponding columns of L(s) are also determined. Specifically, this work focuses on stable and square processes and uses PID controllers; consequently, each ELTF l j (s) does not contain unstable poles, and therefore, closed loop stability is guaranteed if the critical point (−1, 0) is not encircled by the Nyquist diagram of l j (jω). Given a linear margin specification for loop j defined by the line r with parameters j and α j , the PID gains of the j-th column of K(s) must be designed in such a way that every point of the frequency response array l j (jω) to be placed below this line. This condition is defined by inequality (16) at each frequency ω k . This inequality must be parameterized on ρ j , which is the vector of the decision variables, to be used as a constraint into the linear programming problem. Combining (14) with (16) for j = 1 we can rewrite this constraint as follows in (17). Dividing by tangent of α 1 and moving to right the term − (1 − 1 ) provides the desired parameterization, as shown in (18).
This work proposes maximizing the integral gain of the PID controllers to optimize the closed loop performance in terms of load disturbance rejection, since the corresponding integrated error is minimized when integral gains are maximized [34]. Then, for the 2 × 2 case, the optimization problem for loop 1 is formulated as the linear programming problem shown in (19); the same can be defined for loop 2. It is assumed that a previous control matrix K(s) is available to determine l 12 (jω k ) and l 22 (jω k ), and subsequently, χ 1 (jω k ) can be calculated. The inequality constraint in (19-a) corresponds to that in (18). maximize 0 1 0 0 1 0 ρ 1 subject to: The equality constraint (19-b) involves static decoupling for this loop, which implies that non-diagonal elements l ij (0) must be zero at zero frequency. This condition can be parameterized on ρ 1 from (4). In presence of reference changes, static decoupling reduces process interactions at low frequencies to some extent. Furthermore, this constraint considerably improves the convergence of the proposed iterative procedure explained in the next section, since the degrees of freedom are decreased and the search space is narrowed. Other centralized methodologies are mainly focused on obtaining a complete dynamic decoupling performance forcing non-diagonal elements of L(s) to be close to zero at all frequencies. Such imposition usually produces slower disturbance rejection than that achieved with static decoupling [17].
Another issue to note is that the maximization in (19) assumes that the integral gains K I are positive. If this is not the case, the negative signs must be previously removed from ρ 1 and added to the corresponding elements of the χ 1 (jω) array. After optimization, these negative signs must be restored. It is proposed that the signs of the controller elements are derived from the corresponding ones in the inverse of the static gain G(0) of the process. Since the resulting centralized control K(s) achieves static decoupling, K(s) is proportional to G −1 (0), which maintains sign dependency.
The problem in (19) is formulated to obtain PID elements for the control matrix K(s); however, if a PI structure is preferred for some element, its corresponding value of K D in the parameter vector ρ 1 must be set to zero.

C. ITERATIVE PROCEDURE
The optimization problem in (19) is performed independently for each ELTF l j (jω) of the multivariable system to calculate the PID parameter vector ρ j of the j-th column of the control matrix K(s) and to achieve the desired linear margin specification. As it is shown in (14), each ELTF l j (jω) assumes some known elements l ij (jω) of the non j-th columns of the loop transfer matrix L(s), which were calculated from a previous matrix K(s). Therefore, the proposed algorithm needs to define an initial control matrix K 0 (s). Since the optimization problem for each ELTF modifies one column of K(s), the other ELTFs are affected and thus the actual performance obtained by all loops must be verified; if this is not successful, the controller parameters must be updated. Consequently, an iterative tuning procedure is performed. It ends when the design requirements are fulfilled in all loops, or when controller parameters converge within a specified tolerance.
The flow diagram of the proposed algorithm is depicted in Fig. 5. It consists of the initialization of K 0 (s), the calculation of initial open loop transfer functions l 0 ij (jω k ), and the iterative loop through a set of steps performed by each control loop j: calculating each vector χ j (jω k ) at each frequency ω k , tuning of each PID parameter vector ρ j by the optimization in (19), recalculating of all l ij (jω) with the new control matrix K(s), and analysis of a cost index J that checks the design specification achievement. In these steps, the subscript j denotes the loop number from 1 to n, while the superscript m indicates the iteration number.
First, the iteration variable m is set to zero and an initial control matrix K 0 (s) is chosen to calculate the open loop frequency responses l 0 ij (jω k ) before entering the iteration loop. Although the initial K 0 may affect the number of iterations, we have obtained practically the same final control parameters in the analyzed examples when different guesses of K 0 (s) have been used. In this work, the inverse of G(0) is used as the guess for K 0 since it is easily computed, is related with the static decoupling imposed as constraint on the optimization problem, and the procedure seems to reach convergency in fewer iterations when compared with other initial guesses.
Then, for each control loop j and iteration m, the vector χ m j (jω k ) is calculated according to section II.B at each frequency. Next, the corresponding vector ρ m j of PID controller parameters is obtained by the linear programming optimization in (19) to achieve the desired linear margin specification defined by j and α j . Then, the frequency responses l m+1 ij (jω k ) of the new open loop transfer functions are recalculated at every frequency with the new PID parameters, the frequency responses of the updated ELTFs are computed, and the fulfillment of the desired specifications in all loops is evaluated by a cost index J m . If this index is below a user-defined tolerance for three consecutive iterations, the requirements are assumed to be achieved and the design is accepted; otherwise, the iteration variable m is incremented, and the iterative procedure continues calculating the new χ m+1 j (jω k ). When the desired linear margin requirement for a ELTF l j (jω) is evaluated, two possible situations can arise, as shown in Fig. 6: a) all points of l j (jω) are below line r or b) some VOLUME 10, 2022 points are above it. In both cases, the Nyquist diagram must be placed as close as possible to the line r; however, the distance to measure this proximity is calculated differently. In the first case, the shortest distance d min to the line r is selected as the minimum one between the distances of each point l j (jω k ). On the other hand, when some points of l j (jω) are placed above the line r, as depicted in Fig. 6b, the greatest distance d max to the Nyquist plot is selected as the distance. To consider both situations in Fig. 6, the fitting distance D of l j (jω) to line r is defined according to (20). Then, the cost index J for the iterative procedure is proposed as the sum of the fitting distances D j of all the loops, as shown in (21). If this value is below a user-defined tolerance (0.005 per loop), the centralized PID design is accepted.
The main limitation of the proposed methodology lies in not ensuring the convergence. Therefore, the procedure also exits the iterative loop if a maximum number of 10 iterations is achieved, and then, the design associated to the smallest cost index J is selected. Nevertheless, the procedure has worked properly in all tested examples of different dimensions converging after less than ten iterations. Note that there will be always three iterations at least in the proposed procedure.

III. SIMULATION EXAMPLES
Two multivariable processes of different dimensions (2 × 2 and 3 × 3) are simulated to test the proposed methodology. They are classical testbenches from literature presenting multiple time delays and without diagonal dominance, which makes the control design more difficult.
The results are compared with those achieved by other works. The design procedure and simulations are performed with MATLAB R software.

A. EXAMPLE 1: INDUSTRIAL-SCALE POLYMERIZATION REACTOR
This 2 × 2 process [35] is given by the transfer function matrix in (22), where the time scales are in hours. A frequency response array of 1000 logarithmically spaced elements is computed from (22)  The desired linear margins are defined with 1 = 0.64, α 1 = 79 • , 2 = 0.72, and α 2 = 85 • to guarantee lower bounds of 2.75 for A m and 0.63 for M m in the first loop, and 3.5 for A m and 0.72 for M m in the second loop. The proposed procedure obtains the centralized PID control in (23). Different initial guesses of K 0 (s) have been tested: the inverse of the process static gain G R (0), the zero matrix 0 2×2 , the process transfer matrix itself G R (s), the negative allones matrix divided by s, and the multiloop PI controller of Vu [14]. Fig. 7 shows the progression of the PID parameter values for each controller element k ij (s) and the cost index J through the iteration procedure for different K 0 . The inverse of the static gain of the process only needs five iterations for convergence whereas the other guesses require seven iterations. Although these cases differ in the first steps and the number of iterations needed to convergence, all of them achieve the same final control design. Therefore, the initial guess of K 0 has no important effect in the resultant design. However, K 0 = G −1 R (0) usually involves a faster convergence with parameters very close to the final ones after three iterations.
The final resultant Nyquist diagrams of the obtained ELTFs l 1 (jω) and l 2 (jω) are depicted in Fig. 8. Table 1 collects the achieved classical robustness margins and the phase margin crossover frequency ω cp , which is considered close to the loop bandwidth. The proposed method obtains the greatest loop bandwidth. Fig. 9 shows the closed loop system response of the proposed centralized PID control in comparison with other methods, specifically, the centralized PI control of Ghosh [31], the decoupling PID controller of Jin [7], and the multiloop PI controller of Vu [14]. The simulations  are performed with unit step changes at t = 1 h in the first reference and at t = 15 h, in the second one. There is also a load disturbance at t = 30 h that consists of a 0.5 step change in both process inputs. The derivative action of the PID controllers in (23) is filtered by a first order term 1/(|K D /K P |s/N + 1) with N = 20 [34] to achieve realizable implementation and perform the simulation. With this N value, the filter effect can be neglected into the frequency range of interest. Since the proposed methodology does not consider the input sensitivity function, the resultant designs can give large control signals when sudden changes  are applied in the references. In these cases, a first order filter on the references is recommended, such as the filter 1/(0.33s + 1), which is used in this example.
From the simulation, the integrated absolute errors (IAE) and the total variation (TV) of control signals for each loop are calculated according to (24) and collected in Table 1. They are used as performance indices of the control system response for comparison with other methods. The response of the proposed controller shows similar settling times and decoupling performance to set-point changes than other methods that are designed specifically for decoupling. In addition, its load disturbance response is superior since the integral gains of the PID controllers are maximized into the optimization procedure. This better disturbance rejection makes the IAE values of the proposed method smaller than those achieved by the other methods while it has slightly higher TV values. Only the Jin's control obtains a lower IAE value in the first output at the expense of a very conservative response in the second loop with the highest IAE value.

B. EXAMPLE 2: DEPROPANIZER COLUMN
This process [36] is modeled by the 3 × 3 transfer matrix in (25) where time units are in seconds. Its frequency response is approximated by an array of 1000 elements calculated in the frequency range [10 −6 , 10 0 ] rad/s.
In this example, the obtained performance by the proposed methodology is better when the elements of K(s) are forced to be PI controllers instead of PID, and hence K D values are set to zero. The linear margin of each loop is specified by 1 = 0.84 and α 1 = 70 • for the first loop, 2 = 0.66 and α 2 = 68 • in the second one, and 3 = 0.75 and α 3 = 67 • for the third loop. These values are selected according to (15) to obtain similar robustness margins than those achieved by others works for this process. The proposed methodology achieves the centralized PI controller in (26) after five iterations using the inverse of G D (0) as initial guess of K 0 . Fig. 10 shows the progression of the PI parameter values for each controller element k ij (s) and cost index J through the iterative procedure. Although this second example is a higher dimensional system than the 2 × 2 process of example 1, the iteration procedure converges quickly and after two iterations the final design is almost achieved.   Fig. 11 shows the resultant Nyquist diagrams of the three ELTFs l 1 (jω), l 2 (jω) and l 3 (jω), where the fulfillment of the linear margin is verified. In Table 2, the classical robustness margins obtained for each loop are listed.  The proposed centralized PI control is compared with other two centralized designs: the centralized PID control of Garrido [17] and the non-PID control of Wang [36], which are based on decoupling control. The closed loop system responses of the three methodologies are shown in Fig. 12. The simulation is performed with unit step changes in the references at t = 0 s in the first one, at t = 2000 s in the second one, and at t = 4000 s in the third one; there is also a 0.1 step in all process inputs at the same time as load disturbance at t = 6000 s. The proposed control obtains faster responses for reference tracking and load disturbance rejection, particularly in the third loop, where a greater bandwidth is achieved according to the related ω cp value of 0.025 rad/s of Table 2. In addition, the proposed method provides the smallest IAE values with similar or lower TV values than those obtained by the other methods. These values are listed in Table 2.
Only the decoupling performance obtained in y 1 and y 2 with the proposal is a bit worse than that achieved by the other two methods, which are specially intended for decoupling performance. Furthermore, the PI controllers of the proposed design are much simpler than the elements of Wang's control, which are fourth order transfer functions with time delays.

IV. CONCLUSION
A new centralized PID controller design is formulated as an iterative linear optimization problem. The methodology is based on the developed concept of equivalent loop transfer function (ELTF) for centralized control, which allows to design separately the columns of the controller matrix considering the interactions with the other loops through an iterative procedure. The design of each controller column is based on loop shaping for the corresponding ELTF in the Nyquist diagram to satisfy a robustness linear margin. The multivariable PID control design optimizes the load disturbance rejection response of the closed loop system subjected to the robustness linear margin constraint and static decoupling for each loop. The effectiveness of the proposed methodology has been supported by two simulation examples with dimensions 2 × 2 and 3 × 3. Comparisons with other works have illustrated that the proposed method achieves similar or greater performance. The main advantages of the proposed method are the following: • No approximation or model reduction is performed through the iterative procedure, even when the MIMO process has multiple time delays, non-minimum phase zeros or irrational transfer functions.
• The linear margin constraint ensures lower bounds on classical robustness margins in the frequency domain for each ELTF.
• The static decoupling constraint reduces process interactions at low frequencies and allows the method to be applied to non-diagonally dominant systems.
• The method achieves better disturbance rejection responses compared to other centralized control methodologies, such as those based on dynamic decoupling. However, the methodology has some limitations that require further research: • It is restricted to stable processes since the linear margin cannot be applied to unstable systems, where the Nyquist plot must encircle the critical point (−1, 0).
• The convergence of the method cannot be guaranteed, either because of the process characteristics or because the design specifications are not achievable. Nevertheless, the procedure has converged in all tested examples in less than 10 iterations. The static decoupling condition aids considerably to this rapid convergence. The initial guess of the control matrix that is needed to initialize the iteration procedure has little impact on the final design and only affects the number of iterations before convergence.
• The magnitude of the control signals is not considered in the formulation of the method, and consequently, the obtained controller can produce high control values due to sudden changes in the references. In such cases, the references can be filtered with a first order filter.