Design of suboptimal model-matching controllers using squared magnitude function for MIMO linear systems

This paper proposes a novel two-stage method for the design of a suboptimal model-matching controller in an output feedback closed-loop system (OFCLS) using the concept of squared magnitude function (SMF). A streamlined procedure for selection of a reference model, based on a linear quadratic regulator (LQR) with integral action (LQRI) having optimum values for the elements of the weighting matrices and the degree of interaction is proposed. The degrees of the numerator and denominator polynomials of the elements of the OFCLS transfer function matrix (TFM) are obtained from those of the plant and the chosen controller structure. In the first stage of the controller design, taking the LQRI-based closed-loop system (LCLS) as a reference model, the OFCLS is obtained using the approximate model-matching (AMM) technique based on the SMF concept. The approximation method involves a higher-order approximation for stable multiple-input-multiple-output (MIMO) lower-order systems. In the second stage, controller parameters are obtained using the exact model-matching (EMM) method with information about the OFCLS and plant TFMs. The proposed controller design method outperforms the method presented in the literature on integral squared error index. The simulation and experimental results illustrate the effectiveness of the proposed method.


Introduction
The design of controllers for plants with interaction has been one of the challenges faced by researchers in the area of multiple-input-multiple-output (MIMO) control systems. The model-matching technique has been a powerful methodology for controller design of MIMO systems [1][2][3][4][5][6][7][8][9]. Controllers designed using the model-matching technique have been utilized for rejection of various types of disturbances [4]. The concepts of exact model-matching (EMM) of two-dimensional systems using state feedback have been described in [5]. The proposed methodology in [5] can be utilized for the design of two-dimensional filters. A centralized control strategy based on the equivalent transfer function (TF) is presented in [10]. An output feedbackbased controller design methodology for higher-order SISO plants is presented in [11]. The methodology in [11] applies to unstable plant TFs. Model order reduction reduce the complexity of implementation of higher integer-order systems [12,13]. A model-order reduction method has been applied to closed-loop systems with parameter variations [14]. An algorithm for the design of controllers in multivariable systems using the approximate model-matching (AMM) technique has been proposed in [15,16]. EMM of linear time-varying multivariable systems was presented in [17]. A novel EMM methodology for generalized state-space systems via pure proportional state and output feedback was proposed in [18]. EMM has been solved using a generalization of Wolovich's theorem [19]. The modelmatching concept was used earlier for obtaining controllers for tracking and disturbances. In the present work, the AMM technique has been used for obtaining the output feedback closed-loop system (OFCLS) transfer function matrix (TFM), and in a following step, the controller TFM was obtained using the EMM technique.
Trajectory tracking has been achieved successfully using PID controllers in different types of systems [20,21]. Robust adaptive tracking control has been proposed for a nonholonomic mobile manipulator under system parameter uncertainties and external disturbances [22,23]. Trajectory tracking using linear quadratic regulator (LQR) techniques has been successfully applied to different types of systems [24,25]. The weighting matrices of the LQR controller have been determined using a genetic algorithm (GA) [26], particle swarm optimization [27] and pole placement [28]. The methods of determining of weighting matrices described in the literature have not considered the interactions present in MIMO systems. Interaction, as a specification to be preserved in the designed closedloop system, i.e., OFCLS, has also been considered in the present work. The selection of reference model proposed in this work focuses on obtaining the optimal interaction factor for given desired output time-domain specifications together with finding the optimum main diagonal elements of the diagonal weighting matrices to be used in the LQR with an integral action (LQRI)based closed-loop system (LCLS) design procedure. In the present work, the off-diagonal elements of a reference model TFM are determined by considering the optimum interaction present in the controlled system. Squared magnitude continued fraction and factorization techniques are utilized to derive the stabilitypreserved reduced-order approximant [29]. A squared magnitude function (SMF)-based approach has been utilized in the approximation [30] and hardware implementation [31] of single-input-single-output (SISO) fractional-order systems. The retention of stability and minimum-phase characteristics has been achieved in model-order reduction of SISO linear systems using the concept of SMF [32]. The model-order reductions of SISO and MIMO systems have been solved using the method of integration of the stability equation together with the Padé approximation method [33]. The concept of SMF in the field of approximation yields a stable approximant. In the literature, the concept of SMF has been used for obtaining an approximation of SISO fractional-order systems or higher integer-order systems. The present work uses the AMM technique based on the concept of SMF to obtain a higherorder OFCLS representation, subsequently yielding controller parameters.
In the model-matching techniques reported in the literature, optimum interactions present in the OFCLS have not been considered. In [34], the authors developed an algorithm to obtain the desired interaction for a given degree of coupling. In the present work, the design procedure for obtaining a reference model incorporates the time-domain specifications and optimal interaction utilizing a simultaneous GA-based optimization procedure.
The rest of this paper is organized as follows. Section 2 presents the problem statement. The proposed method for selection of the reference model is described in Section 3. Section 4 details the proposed AMM algorithm based on the SMF concept. Section 5 describes the procedure to determine the controller parameters based on the EMM technique. Simulation and experimental results are provided in Sections 6 and 7 respectively. Finally, the concluding remarks and future scope are included in Section 8.

Statement of the problem
The objective of the present work is to design a controller G c (s) for a MIMO plant G P (s) using the modelmatching technique. The controller design aims to make the response y(t) of the OFCLS G CL (s) as close as possible to the response y d (t) of the reference model G d (s). G d (s) is chosen in such a way that it reflects the desired closed-loop system characteristics, such as the desired time-domain specifications and the optimal interaction to be present in the G CL (s). Model-matching controller design in the literature usually deals only with one of the two types of model-matching, namely EMM and AMM. The present work deals with both of these model-matching techniques for the design of G c (s), thereby combining the advantages of both types into one algorithm. The important aspect of the controller design in the design and implementation is that the controller structure need not be fixed. The user can choose any controller structure to achieve the objectives. In addition, the solution for the controller parameters need not be unique. Certain controller parameters can be chosen, and the remaining controller parameters can be obtained according to the developed algorithm, resulting in a reduction in hardware compulsion. The present work leads to the retention of the stability property of the G d (s) in the G CL (s), thereby yielding a stabilizing G c (s) in the second stage of the algorithm. The block diagrams of the reference model (LCLS model/user-defined model) and the OFCLS are shown in Figure 1.
The plant can be represented in the s-domain as where G . . , q represents the elements of the plant TFM relating the kth input to the ith output. The numerator and denominator polynomials of each element of the plant TFM are represented by n (i,k) P (s) and d P (s) respectively. The problem of finding G c (s) such that G CL (s) yields the same response as that of G d (s) can be stated as follows. Find G c (s) for the given MIMO plant from the information about the degrees of the numerator and denominator polynomials of elements of G c (s) with the objective of obtaining the optimal interaction and specified time-domain specification in G CL (s). In the proposed method for controller design, the first stage consists of an approximation of G d (s) to a higherorder TFM G CL (s). The degrees of the numerator and denominator polynomials of each element of G CL (s) are obtained from those of the plant and the chosen controller structure (see Appendix A). In the second stage of the proposed method, G c (s) is obtained by using the concept of EMM with the knowledge of G P (s) and G CL (s).

Reference model selection
This section proposes a streamlined procedure for the selection of G d (s) which reflects the optimum degree of interaction and embodies the desired time-domain specifications.

Conventional model selection approach
Methods based on the preservation of the pole-zero excess of the plant [2] and specified time-domain specifications [34] have previously been utilized for G d (s) selection. The approach presented in [34] embodies a specified interaction in the off-diagonal elements of the G d (s) TFM. The limitation of the G d (s) selection approaches presented in the literature is that they do not include information about the plant dynamics.

Proposed model selection procedure
The proposed method uses the LCLS model as G d (s) for the design of G c (s). The method uses an initial model which is utilized to compute the LCLS model. The initial model is framed with the desired time-domain specifications. The interaction parameters λ (i,k) of the initial model are retained as tuning parameters. The LCLS model is framed by tuning the elements of the weighting matrices Q and R in the optimal state feedback design procedure. The interaction thus developed in the LCLS model is mapped onto the initial model by tuning the parameters λ (i,k) . The desired time-domain specifications of the initial model are used to tune the elements of the weighting matrices of the LCLS model. Thus, the LCLS model reflects the desired time-domain specifications and the optimal interaction to be present in G CL (s).

Initial model
The desired transient response specifications include natural frequency, settling time, damping factor, etc. The initial model dynamics can be expressed in the form of TFM as Here, each G (i,k) m (s), i = 1, 2, . . . , p; k = 1, 2, . . . , q represents an element of the initial model TFM relating the ith output to the kth input. The Laplace transform of the desired responses is denoted by Y (i) m (s), with steady-state values denoted by S i . The desired profiles of Y (i) m (s) can be taken in the form of the response of a standard second-order lowpass filter for R i (s), taken as a step signal of amplitude S i , i.e., where ω n is the natural frequency and ζ is the damping ratio, both of which can be determined from the desired time-domain specifications such as settling time, overshoot, rise time, etc. The (i,k)th off-diagonal elements of the initial model TFM are designed to have a pole at s = −λ (i,k) and zero at s = ∞, which is utilized along with Equation (3) in Equation (2) to yield the expressions for diagonal elements of the initial model: where u = 1, 2, . . . , min(p,q); k = 1, 2, . . . , q. For a given desired time-domain specification, the dynamics of G m (s) will be a function of the parameter λ (i,k) . A given value of λ (i,k) quantifies the interaction that can be present in the controlled dynamics. Higher values of λ (i,k) imply a system with lesser interactions and vice versa. In the present work, a method of obtaining the optimized value of λ (i,k) along with optimized values of the elements of the weighting matrices, is presented in subsection 3.2.3.

LCLS model
To eliminate the steady-state error, the integrating controller has full state feedback incorporated. The difference between the reference input and the output is integrated and the outputs of the integrators are considered as additional states. The combined structure is as given in [35].
In the combined structure, x ∈ R n×1 , u ∈ R q×1 , y, x i , r ∈ R p×1 are the state, input, output, additional state and reference vectors respectively. The numbers of inputs, outputs and states present in the plant dynamic model are denoted by q, p and n respectively. The state-space representation of the plant G P (s) shown in Equation (1) can be given aṡ where the constant matrices A ∈ R n×n , B ∈ R n×q , C ∈ R p×n , K ∈ R q×n , and K i ∈ R q×p are the state, input, output, state feedback gain and integral gain matrices respectively. Elimination of steady-state offset is achieved by incorporation of the integrator into the controlled system [35]. The overall TFM representing the LCLS model can be obtained by where the optimal state feedback gain matrix is given by . K op is obtained by solving the algebraic Riccati equation for a given weighting matrices Q and R [25,26]. In the present work, the weighting matrices are obtained using the optimization procedure explained in the following subsection.

Weighting matrices and degree of interaction
The set of optimal weighting matrices and the degree of interaction have been determined by minimizing the difference in step response of G m (s) and G d (s). In the present work, the weighting matrices Q and R of the LQRI are assumed to be diagonal matrices. Hence, the aim of the optimization process is defined as follows.
Find the main diagonal elements of diagonal weighting matrices of LQRI and the off-diagonal pole location of the initial model TFM to minimize subject to the constraints that the weighting matrices Q and R should be positive semi-definite, and positive definite matrices, respectively.
where the step responses of TFs relating the ith output to the kth input of G m (s) and G d (s) are represented by y (i,k) m (·) and y (i,k) d (·) respectively. The sample time and final time for the evaluation of the step response are t and τ t respectively. The superscript T denotes the transpose of the matrix.
The solution of the optimization problem leads to an optimized value of λ (i,k) and optimized values of the elements of the weighting matrices. The selection of the LCLS model over the initial model for G d (s) as an input to the first stage of the proposed modelmatching-based controller design procedure guarantees two major advantages. Firstly, the LCLS model is formulated by considering the dynamics of the plant. Secondly, the degrees of the numerator and denominator polynomials of elements of the LCLS model TFM will be closer to those of the OFCLS TFM in comparison to those of initial model TFM. The overall procedure for obtaining the optimal LCLS model from the information in the initial model can be illustrated in the form of a flow diagram as shown in Figure 2

AMM technique based on SMF concept
This section presents the development of the proposed approximation algorithm for MIMO systems based on the concept of SMF. As part of the first stage of the proposed algorithm, an LCLS model TFM is approximated to the OFCLS TFM using SMF. If the degrees of the numerator and denominator polynomials of elements of the LCLS model TFM are same as those of the OFCLS TFM, the G c (s) parameters may be obtained directly using EMM by taking the LCLS model itself as the OFCLS. In all other cases, the SMF-based approximation plays a role in the determination of OFCLS TFM.

Formulation of non-homogeneous simultaneous equations
This section describes the procedure for constructing non-homogeneous simultaneous equations, which are solved with the objective of finding parameters of SMF of OFCLS. Let G d (s), obtained in the form of TFM via Equation (6), be expressed as . . , q represents the elements of the G d (s) TFM relating the kth input to the ith output. Elements of the TFM of G d (s) can be represented as and b β , β = 0, 1, . . . , n 1 are the numerator and denominator polynomial coefficients of the (i,k)th element of the LCLS model TFM respectively, having m (i,k) where G represents the elements of OFCLS TFM relating the kth input to the ith output. The numerator and denominator polynomials of each element of the OFCLS TFM are represented by n (i,k) CL (s) and d CL (s), respectively. Each element of the OFCLS TFM can be represented as and c μ , μ = 0, 1, . . . , n − 1 are the numerator and denominator polynomial coefficients of the (i,k)th element of the OFCLS TFM respectively, having m (i,k) ≤ n. Each element in the SMF of the OFCLS TFM, P (i,k) CL (s 2 ), can be represented as The right-hand side of Equation (12) can be represented as a ratio of two even polynomials. Hence, The numerator and denominator polynomial coefficients of P and c μ as given in Equation (11). The OFCLS TFM is designed with the objective that the frequency response of its SMF matches that of the LCLS model at certain frequency points, s = jω h , within a desired frequency range. i.e., where P (i,k) d (s 2 ) is the SMF of G d (s) and ω h , h = 1, 2, . . . , n x , are n x frequency points along the positive imaginary axis in the s-plane. The value of n x is computed using a ceiling function, as where n nd is the number of the unknown numerator and denominator coefficients in the SMF of the OFCLS TFM. The numerator polynomials of the SMF of the OFCLS in Equation (13) can be denoted in matrix form as where S Similarly the denominator polynomial in Equation (13) can be denoted in matrix form as where S C = [s 2(n−1) . . . s 2 1], and X C Then, combining Equations (14), (16) and (18), the following matrix equation can be formed: where h = 1, 2, . . . , n x . In general, for each frequency point ω h , Equation (20) leads to the formation of pq equations. The pq set of matrix equations can be cascaded in the form where For a 2 × 2 MIMO system, the matrices A h and B h can be shown as in Equations (23) and (24) respectively.
where O 1 ×(m (i,k) +1) represents the null matrix of dimension 1 × (m (i,k) + 1), i = 1, 2, k = 1, 2 and The matrices A h and B h obtained at each frequency point are cascaded to form the resultant A and B matrices respectively. The resulting set of non-homogeneous simultaneous linear equations are obtained in the form The least squares solution of Equation (25) yields the values for the unknown parameters in vector X as given in Equations (17), (19) and (22), thereby obtaining the SMF of OFCLS.

Closed-Loop system steady-state matching
The present work considers the steady-state matching of G d (s) with G CL (s), while the input signal is a unit step using the relation Let f given in Equation (14) at the origin of the s-plane. Then, equating the value of SMF of the OFCLS at the origin of the s-plane to that of the LCLS using Equations (13) and (14) respectively, yields As a result of Equation (27), the number of unknown coefficients n nd to be determined is reduced by the product pq for a p × q MIMO system. Hence the total number of frequency points to be optimized during the approximation procedure is also reduced. In the present work, Equation (27) is incorporated into Equation (25) to ensure steady-state matching between G d (s) and G CL (s), while the input signal is a unit step signal.

Factorization and selection of roots of the numerator and denominator polynomials of SMF
The numerator and denominator polynomials of the SMF of OFCLS P (i,k) CL (s 2 ) will be distributed symmetrically with respect to the real and imaginary axes of the s-plane. The poles and zeros of the SMF of the OFCLS are determined via the factorization of polynomials n (i,k) CL (s 2 ) and d CL (s 2 ). The present work incorporates the possibility of complex poles and complex zeros of P CL (s 2 ). The real and complex conjugate pair roots of P CL (s 2 ) present in the LHP are retained in the OFCLS using Equation (28), to form the real coefficients in the numerator and/or denominator polynomials of G (i,k) CL (s).
where m are on the imaginary axis of the s-plane. The present work tackles the problem of having purely complex roots of the numerator and/or denominator polynomials of P CL (s 2 ) with the help of an optimization framework. The objective of optimization is to minimize an objective function subject to the constraint that the set of chosen frequency points avoids the case in which any poles/zeros of P CL (s 2 ) become purely imaginary.

Selection of optimal frequency points
The problem of choosing optimal frequency points is formulated in an optimization framework where the   fitness function is constructed by calculating the sum of squared differences in the corresponding real and imaginary parts of the frequency responses between the LCLS model and the OFCLS, within the desired bandwidth. Hence, the aim of the optimization process is defined as follows.
Find the frequency points, ω h , h = 1, 2, 3, . . . n x to minimize 29) where ω r , r = 1, 2, . . . , M are equally spaced frequency points within the desired frequency range of interest. Re{·} and Im{·} denote the real part and imaginary part respectively, subject to the constraints that the roots of the numerator and/or denominator polynomials of the SMF of the OFCLS should not lie on the imaginary axis of the s-plane. The optimization problem presented in Equations (7) and (29) was solved using GA.

Determination of controller parameters based on EMM technique
This section presents the procedure for determining controller parameters based on the EMM method with knowledge of the OFCLS and plant TFMs. Let G c (s) be represented in TF domain as where G i = 1, 2, . . . , p; k = 1, 2, . . . , q represents elements of G c (s) relating the kth controller input to the ith controller output. The numerator and denominator polynomials of each element of G c (s) are represented by n (i,k) c (s) and d c (s) respectively. The OFCLS TFM shown in Figure 1 can be expressed as a function of G p (s) and G c (s) as

Example 1
The proposed suboptimal controller design method is illustrated for the design of a model-matching controller with the objective of tracking the desired circular trajectory of a wheeled mobile robot (WMR), namely, QBot 2. Details of the modelling of the WMR dynamics are as presented in [34].

Determination of parameters of the initial model and LCLS model for WMR
This section presents the development of optimal parameters of the initial and LCLS models for the WMR. The settling time and ζ of the wheel velocity profiles of the WMR are taken as 3 s and 1.1 respectively. In the present work, S 1 and the desired radius of the circular trajectory are set at 0.2 m/s and 0.3525 m respectively. The value of S 2 is then calculated as 0.1 m/s [34]. Directions for implementing the optimization technique in the determination of the LCLS model using the global optimization toolbox of MATLAB are presented in Appendix B.
The comparisons of the right and left wheel velocity profiles due to the step of 0.2 m/s and the step of 0.1 m/s respectively of an initial model to those of the tuned LCLS model, are presented in Figure 3(a,b). The  comparison of the resultant trajectories of the WMR is presented in Figure 3(c). Figure 4 shows the comparison of the step responses of the optimal LCLS model with those of the initial model. From the Figure 4 it can be seen that the desired time-domain specification has been embodied in the LCLS model with weighting matrices Q and R chosen optimally. The response of the off-diagonal elements of the LCLS model TFM is close to that of the corresponding elements of initial model TFM with a maximum error of the order of 10 −3 m/s. The present work evaluates the objective function in Equation (7) during a time interval T 2 , which is the time taken by the WMR to complete two revolutions along the desired circular trajectory. The parameter τ can be calculated as τ = T 2 /( t). In the present work, t is taken as 10 −2 s. The time taken for two revolutions of the WMR along the desired trajectory is calculated as 29.531 s [36].
The tuning of the weighting matrices presented in [25,26] is performed with the objective of mapping the desired time-domain specifications as a certain combination of weights of states and inputs. In the present work, the tuning of the weighting matrices considers the mapping of the interaction present in the initial model onto the LCLS model, in addition to the desired time-domain specifications.

Determination of OFCLS TFM
The present work illustrates the proposed algorithm for different G c (s) structures as given in Table 1. The design of a practical PID using the proposed method is carried out by incorporating a low-pass filter along with the differentiator [37][38][39].
The LCLS model is approximated to a TFM having the degrees of the numerator and the denominator polynomials as decided by the chosen controller structure, as shown in Table 1. In the present work, the desired frequency band of approximation and the value of M are taken as (10 −3 −10 2 ) rad/s and 700 respectively. Optimal frequency points obtained in the case of each chosen controller structure are given in Appendix C. Equation (32) is utilized in the EMM technique to obtain the unknown controller parameters (see Appendix D).

Comparison of frequency responses
The performance index J 2 in the first stage was minimized with the objective of matching the frequency response of the OFCLS to that of the LCLS model. The set of figures given in the first column of Figure 5 shows the comparison of the frequency response of the OFCLS with that of the desired LCLS model for PID-based controller structures, whereas the set of figures in the second column shows the comparison with lag-based controller structures.
From Figure 5, it can be concluded that the frequency responses of G d (s) and G CL (s) match well. Figure 6(a,b) show the comparison of the trajectory of the WMR obtained by OFCLS TFM with that obtained with the desired LCLS model TFM, for PID-based G c (s) structures and lag-based G c (s) structures respectively. Figure 6(c,d) show the variation of the x-coordinates of the WMR for PID and lag-based G c (s) structures respectively. Figure 6(e,f) show the variation of the ycoordinates of the WMR for PID and lag-based G c (s) structures respectively. Figure 6 shows that the actual trajectory of the WMR using the designed controllers satisfactorily matches the desired trajectory. The computational time taken for each stage of the algorithm is shown in Table 2. The computational times shown are for MATLAB R2015b software running on a 3.3 GHz Intel(R) Core(TM) i5-4590 processor with 8 GB RAM.

Example 2
The proposed method has been utilized to design a decentralized G c (s) for a 10th order MIMO plant taken from the literature [2]. The G d (s) in the first stage of the algorithm for this example has been selected based on pole-zero excess as described in [2]. With this G d (s), the algorithm is started using the AMM presented in Section 4. The step response of the G d (s) is compared with that of G CL (s), obtained using the proposed method and the method presented in the literature [2], as shown in Figure 7. Comparisons of the minimum value of the objective function (Equation (29)) and maximum absolute deviation obtained using the method presented in the literature [2] with the values obtained using the proposed method, are presented in Table 3. Table 3 shows that the proposed method of controller design is better in comparison to the method reported in [2] on integral squared error index. The G c (s) parameters obtained using the proposed method are presented in Appendix 5.

Experimental results: Example 1
The operation of the WMR is based on the communication between the host computer and the target, QBot 2 [40]. The G c (s) having the chosen structure and designed using the proposed algorithm is implemented on the host computer using the MATLABbased Simulink software. The real-time control software, QUARC downloads real-time code generated from the host computer to the embedded computer mounted on the QBot 2 platform [40]. Figure 8(a,b) show the comparisons of the trajectory of the WMR obtained by OFCLS TFM with that obtained using the desired LCLS model TFM for PID-based G c (s) structures and lag-based G c (s) structures respectively, obtained in an experimental set-up. Figure 8(c,d) show the variation of x-coordinates of WMR for PID-based G c (s) structures and lag-based G c (s) structures respectively, obtained in an experimental set-up. Figure 8(e,f) show the variation of the y-coordinates of WMR for PID-based G c (s) structures and lag-based G c (s) structures respectively. Snapshots of the experimental output are presented in Figure 9(a-e). The blue circular trajectory is the desired trajectory, and the centre of the axle of the WMR tracks this trajectory. The deviations present in the experimental results are greater than those of the simulation results; this may be due to modelling errors incorporated during formulation of the mathematical representation of the plant dynamics.

Conclusion
In this paper, a novel methodology for the design of suboptimal model-matching controllers for MIMO systems based on the SMF concept is presented. The conventional formulation of the reference model depends on the experience and prior knowledge of the designer about the requirement and performance of the plant. The advantage of the proposed reference model selection is that it involves the dynamics of the plant. The design procedure during the selection of the reference model encounters two problems. Firstly, interactions to be present in the designed closed-loop system are difficult to quantify as a specification. Secondly, in the LQRI design procedure, it is difficult to associate the desired time-or frequency-domain specifications with the elements of weighting matrices Q and R. The proposed reference model selection procedure addresses these two problems by formulating the problem as an optimization problem. The developed controller design algorithm is also compatible with user-defined reference models as illustrated in the second example.
As a part of the procedure in the first stage, a higherorder approximation method for MIMO lower-order systems based on the SMF concept is proposed. The LCLS model requires information regarding all the states for the purposes of implementation. The performance of the designed higher-order OFCLS mimics that of the lower-order LCLS model. Thus, the designed OFCLS can be termed suboptimal. The OFCLS is designed to have same denominator polynomial in all the elements of its TFM which results in its minimum realization. The developed algorithm also incorporates the complex poles and/or zeros in the SMF of an approximant. The SMF-based approximation ensures the steady-state matching of the reference model and OFCLS. The approximation procedure also ensures the stability of the OFCLS. Analysis of the closed-loop system can be performed at the end of the first stage itself.
In the second stage, the controller parameters are obtained using the EMM method with the information of the OFCLS and plant TFMs. In the EMM method, some of the unknown controller parameters can be assumed, taking into consideration the availability of hardware components during implementation. The rest of the controller parameters can be obtained by the solution of the resulting simultaneous equations. In this work, the filter coefficient present in the practical PID is kept constant at a value of 100, and the rest of the parameters are obtained as described. The proposed method for the controller design is illustrated by taking two examples from the literature. Work has also been carried out on utilizing the proposed algorithm for the design of fractional-order controllers for MIMO systems. This work will be reported elsewhere.

Disclosure statement
No potential conflict of interest was reported by the author(s).  Pc + m (2,1) Pc , n Pc + m (2,2) Pc ), and Deg

AB.1. Settings of GA
The weighting matrices Q ∈ R (n+p)×(n+p) and R ∈ R q×q of LQRI were taken as the diagonal matrices. The dynamic model of the plant [34] had n ( = 4) states and together with the p ( = 2) additional states of LQRI, the total states present in the LQRI becomes n + p ( = 6). Hence the number of parameters to be optimized for the weighting matrices Q and R were obtained as 8. Together with the pole location of off-diagonal elements of the initial model, the total number of parameters to be optimized nvars becomes nine. The set of optimal parameters ν (1) op has been chosen by utilizing MATLAB function "ga". The optimal parameters obtained in this paper were obtained by selecting the lower and upper bound for the weighting matrices Q and R as zero and infinity respectively. The lower and upper bounds for the degree of interaction λ were set as 100 and 4500 respectively.

AB.2. Determination of the initial model and the LCLS model from ν (1) op
The optimized values of the initial model and the LCLS model parameters obtained by using GA are λ = 3167. 6

Appendix C. Direction for Obtaining OFCLS: Example 1
The lower and upper bounds for the frequency points were determined as presented in [31]. The optimal frequency points ν (2) op obtained during the approximation procedure for different chosen controller structures are presented in Table  A1.
The determination of ν (2) op is followed by finding the solution of Equation (25), to obtain the coefficients of the SMF of OFCLS TFM. Now the stability preserving root selection strategy as mentioned in subsection 5.3 can be utilized to obtain the OFCLS TFM.