Computing frequency response of non-parametric uncertainty model of MIMO systems using υ -gap metric optimization

This paper presents the computation of the non-parametric uncertainty model for multi input multi output (MIMO) systems, which is described by normalized coprime factors (NCF) using the frequency response data of the system. This computation is accomplished by minimizing a υ -gap metric criterion. For this purpose, the problem is formulated to a convex optimization context, such that a semideﬁnite programming (SDP) can be imple-mented. Minimization constraints and the normality constraints of coprime factors are converted to linear matrix inequalities (LMI). Thus, by convex optimization algorithms, the semideﬁnite programming will be optimized. The proposed method can also be used for non-square multi input multi output systems in a conservative assumption. So, through the ﬁrst process of optimization, the frequency responses of the normalized coprime factors are derived. Finally, to evaluate the performance of the proposed method in the computation of the normalized coprime factors of a system, the simulated results of this method are compared with those obtained by the other methods for two types of systems.


INTRODUCTION
For years, control-oriented system identification techniques have been presented to deliver the proper models including a nominal model and a form of uncertainty [1]. The behaviour of perturbed systems may be different when they are measured by an operator norm under feedback control. The models used in many robust control techniques, like loop shaping, are required to be in the form of normalized right or left coprime factors and direct identification or computation of these factors is very important and applicable. Based on frequency response data direct identification of these factors is possible by computation of the frequency response of coprime factors. On the other hand, the gap metric criterion is used to show that the problem of robust optimization is equivalent to robust optimization for normalized coprime factors (NCF) perturbations [2]. It is well known that a relatively small uncertainty in lightly damped poles and zeros can result in a large distance measured in the ν-gap metric, leading to conservative robust stability and performance guarantees [3]. However, while many uncertainty models with additive and multiplicative forms have been presented, the NCFs model has more appreciable properties in rep- resenting the system dynamics, especially for an unstable system [4]. Since there is a clear relationship between the υ-gap metric and NCF models, this model is more attractive for perturbed system identifications. However, compared to the identification of an additive perturbed transfer matrix model, more efforts are required to identify the NCFs model from experimental data. This is because the NCFs model is a purely mathematical representation of plant dynamics, and their time or frequency domain response cannot be directly measured [5].
In some recently researches, some methods have been introduced for the identification of the NCFs model, especially in a closed-loop for unstable systems, where it is required to have a controller in the loop [6], so, model identification and controller design are accomplished simultaneously. In such research and similar works, knowledge of the initial controller is needed for identification, and a certain form of NCF is identified with parameterization. The main goal of these researches, such as [4,[6][7][8][9][10], is to find an optimized controller that can minimize the performance index. Moreover, some presented methods are only applicable to SISO systems. Some other methods, like ones that use the orthogonal basis functions, cannot be applied to all systems, especially unstable systems [11]. The methods in [5,12] and, which are close to our work, are not as straightforward as the method presented in our paper. In these two methods it is required to use a parametric model (matrix) to ensure convergence of the algorithms. In addition to the stability of the nominal model, normality and coprimeness are needed to identify NCFs. To overcome these difficulties, it is suggested in [13] that in order to identify an NCFs model, it may be appropriate to, at first, estimate the time or frequency domain response of NCFs from measurable physical variables. When this estimate and its error characteristics are available, the estimation of a plant's NCFs and validation of a priori information about them are possible. In [13], NCFs are identified with constrained curve fitting. But, in practice there is not a mathematically appreciable relation between NCFs and physical measurable variables [5].
However, based on our best knowledge, until now, the presented methods for identification of MIMO systems based on the general form of the perturbed NCFs model, are challenging. Generally, when an indirect approach is adopted in nonparametric estimation, the estimate is possibly with an infinite variance; when a direct approach is employed in deriving a parametric model, a noise model is required that can describe the actual disturbance characteristics etc. [5]. For example, [14] presented robust controller design and identification of the coprime factorization model for MIMO systems. In order to convexify the problem, they used Taylor expansion at the initial point, considering a closed-loop stable configuration for a given initial controller. A newer one [15] developed an FRD control design for stable or stabilized MIMO for any kind of control structure based on non-smooth optimization.
However, they (and other authors) identified coprime factors (in a parametric method) through control design, but these methods are out of our scope and we did not add them to our references and literature surveys.
In this paper, we will show that based on the υ-gap metric and frequency responses of a system, we can introduce a semidefinite programming (SDP) whose output variables are right or left NCFs and the resulted models have the least difference with frequency responses according to υ-gap metric. So, in this new method, by a standard form of SDP, frequency responses of NCFs are computed. After computing the frequency response of the NCFs with presented methods (as in [13]), the NCFs model will be identified, showing that the derived nominal model has the least difference with a real model.
After preliminaries (Section 1.2), in Section 2, we will discuss mathematical problem formulation and extract a formulation based on frequency responses for the NCFs model and υ-gap metric. In Section 3, the proposed approach for the computation of NCFs in the frequency domain is provided and the results are summarized in a main theorem. Moreover, at the end of this section, the steps for the problem solving, which is a common SDP problem, are provided in a proposed algorithm. In the Section 4, we discuss the simulation of two examples and compare the results of the proposed algorithm with the ones in the two references by some related graphs. Finally, in Section 5, we conclude the paper. In the annex, some relations and proofs which are needed for discussion, are delivered.

Notations
In this paper, capital and bold letters denote matrices, while small and bold letters denote vectors. For scalars, regular alphabets are used. We use * for the conjugate transpose of matrix. If we use ∼ in front of a matrix in the s space, it means A ∼ (s) = A T (−s). Also, for a matrix A − * = (A −1 ) * and |A| represents the determinant of A.
To show that a matrix A is positive semi-definite or positive definite, we use A ≽ 0 or A ≻ 0 respectively; for negative semidefinite or negative definite cases, we use ≼ or ≺ respectively. If we have A ≽ B it means A − B ≽ 0.

Preliminaries
The transfer function matrix of any MIMO system with n inputs and m outputs is an m × n matrix, which can be written in the NCFs form, as follows.
Here, M mm (s), N mn (s) ∈  ∞ are left coprime factors(LCFs) and P mn (s), Q nn (s) ∈  ∞ are right coprime factors (RCFs). If these matrices (factors) are normalized, they will satisfy P * For the perturbed MIMO system, the normalized right coprime factorization is defined as where P * (s)P(s) + Q * (s)Q(s) = I, and W P (s) and W Q (s) are weighting functions matrices and [ Δ P Δ Q ] is the uncertainty matrix. Also, the normalized left coprime factorization can be defined similarly for perturbed MIMO systems that also is used in this paper.
On the other hand, in the robust control context, υ-gap metric is used for the determination of the maximum distance between two transfer functions; it is a useful tool in uncertainty analysis and feedback system's robustness. To calculate the υ-gap metric, assume that G mn is the transfer function matrix of an MIMO system described by RCFs and LCFs as G(s) = P(s)Q −1 (s) = M −1 (s)N(s) which can be represented in the frequency domain. Also, assume G i is the transfer function matrix of another system described by RCFs and LCFs as The υ-metric criterion between these two systems, is defined as [16] (G, In this definition, wno is the abbreviation representing the number of counter clockwise encirclements around origin by | | asevaluated on the Nyquist contour [16].
In this paper, three following theorems are used.

Theorem 1. [17] Given the complex matrices A, B and C of compatible dimensions, there exists a solution Δ to the linear matrix equation
if and only if A − BC −1 B * ≽ 0.

Theorem 3.
[17] LetF = F * , L, R and D be the complex matrices of appropriate sizes, with D ≤ 1. Let Δ be a complex perturbation that has the block diagonal structureΔ = diag( Δ 1 , … Δ N ). Partition L, R and D conformally with Δ as: Then, a sufficient condition for

PROBLEM FORMULATION
In this section, the problem of computing the general form of the NCFs model for the MIMO perturbed system in the frequency domain is formulated, which is expressed in a problem of semi-definite programming (SDP) with a υ-gap metric index. This problem can be solved by convex optimization algorithms.
Assuming the frequency responses of the system G are available, we are going to obtain the frequency response of the NCFs of system. For this, let us have the frequency responses of a system in some different frequencies ( 1 , … , N ), which are taken at different operating conditions or of different samples of a product. Consequently, {G i ( j )} is a set of frequency responses of the system that are measured in the frequency j and i = 1 … K conditions. So, in different operating conditions or for several system samples, there are different system models; therefore, it is needed to have a nominal and an uncertainty model whose behaviour has the least difference with available responses (G Δ = {G 1 , ⋯ ,G K }), according to some indexes. So, we will define this difference with υ-gap metric. In mathematical expressions, we want to determine the normalized right coprime factors (NRCFs) in a general form as follows.
It is needed to compute the behaviour of factors P and Q and their weighting matrices for all frequencies ( 1 , … , N ).So, in the worst condition, the distance between the nominal model (G) and G Δ is minimized based on the υ-gap metric. The worst condition of the distance between the nominal model and G Δ set, according to the υ-gap metric, is defined as: Hence, the nominal model which has the least J(G) is the best [18]. According to the υ-gap metric, minimization of J(G) means It is clear that if in each frequency there is max , the above minimization constraint will be satisfied, which means that in each frequency, we have the following condition [19] max So, if this minimization is done in each frequency, the amount of difference orJ(G) index will be less than max = max 1 ,…, N . Consequently, by solving the above optimization problem, not only the frequency response of the nominal model of G, but also the weighting functions matrices (W P (s) and W Q (s)) and the υ-gap metric of (G,G Δ ) = √ will be obtained. Since we want to compute the frequency response of NCFs (P and Q), the formulation will be rewritten according to these factors. This formulation is done for the left NCFs in the following.
Using the definition Equation (4) for υ-gap metric and the following definitions where

MAIN RESULTS
According to the formulation (10), during the minimization problem in each frequency, will be minimized and X 0 ( j ) and W( j ) will be computed. To compute the frequency response of a model which has the minimum , the problem is converted to a minimization problem with matrix inequality constraints by presenting the following lemmas.

Lemma 1. If the following matrix inequality is satisfied
Proof. By using the Theorem 2, Equation (11) can be converted to the following relation. For simplicity, So the minimization constraint J(G) in Equation (10) is converted to the above constraint, which of course, is not linear. At each frequency, the constraint Equation (13) can be changed to the following linear matrix inequality (LMI) So, according to the Theorem 3, a sufficient condition for having Equation (15) is that there exists a matrix = diag( 1 I, … , n+m I), partitioned conformally with Δ, with i real scalars, such that the inequalities in Equation (13) are approved. □ However, the matrix X 0 should have a special property that is made by right NCFs and satisfies X * 0 X 0 = I, which is defined for a unitary matrix. To define some LMI constraints to guarantee the unitary condition of X 0 , the following lemma is presented.

Lemma 2. For semi-definite matrices A and C, and the invertible B with appropriate sizes
Proof. According to the Theorem 2, if is positive semidefinite, there is A ≽ X * B −1 X. On the other hand, because of C ≽ A ≽ 0, there is C ≽ X * B −1 X; because the condition in the Theorem 2 is sufficient and necessary, then The normality constraint of the coprime factors is equivalent to a unitary constraint of X 0 , which is in Equation (10). In many optimization algorithms, there is an iteration loop and in each step, the previous values of variables are used in the new step. So, in our SDP optimization loop, in each step, X 0 will be computed and used in the next step. We want to minimize the cost function and have the new value of X 0 converge to a unitary matrix simultaneously. For this purpose, the following two lemmas will be proved for the guarantee of convergence X 0 to a unitary matrix in an iterative loop.
Proof. If X 0 is a unitary matrix then X * 0 X 0 = I which is not a convex constraint. According to the Theorem 1, the constraint guarantees that there is a contractive matrix Δ, such that On the other hand, for any matrixX we have So, So, if we addX * X 0 + X * 0X −X * X ≽ I then it satisfies that X * 0 X 0 ≽I. On the other hand, for any matrices A and B, if ≽ 0 [20]. So we can rewrite as If Hence, according to the Lemma 2, it ensures the Equation (18). Finally, this constraint ensures Equation (19). So, in an optimization loop, we try to have the contractive matrix Δ reach a unit matrix and as a result, the matrix X 0 will be unitary.
Assume that 0 < ≪ 1; so, the condition of results in the following BMI as According to the Theorem 2, the Equation (24) can be represented as Equation (17), which is an LMI constraint. If the Equation (16) is added to the constraints of the problem and the matrixX is the matrix value of X 0 of the previous step of the optimization loop, we have at first Also, sincē(Δ) ≤ 1, then 1 ≤̄(X 0 ). So, this constraint with the constraint in Equation (17) ensures that Hence, in an iterative loop, the matrix X 0 converges to a unitary matrix with an accuracy of 0 < ≪ 1. However, because of 0 < ≪ 1, the matrix X 0 can be approximated as a unitary matrix. □ Also, to guarantee that the matrix X 0 converges to a unitary matrix and the contractive matrix Δ to a unit matrix, another condition is needed.

Lemma 4. The following LMI
guarantees that 1 ≤̄(X 0 ) ≤̄(X); in an iterative loop, and the matrix Δ in Equation (19) converges to a unit matrix.
Proof. If we have the conditionX * X ≽ X * 0 X 0 , it guarantees that 1 ≤̄(X 0 ) ≤̄(X); in an iterative loop, thē(X 0 ) will be small and smaller and finally, it tends to be 1. In this situation, the matrix X 0 also tends to a unitary matrix. According to the Theorem 2, the conditionX * X ≽ X * 0 X 0 can be written as the LMI in Equation (27). This LMI will be true if and only ifX * X ≽ X * 0 X 0 and consequently,̄(X 0 ) ≤̄(X). Since, for any non-singular matrices A or B, with proper dimensions, there arē Therefore, with the assumption of Equation (17), there will be So, according to Equation (27), in each iteration there is (X 0 ) ≤̄(X); according to Equations (28) and (29), the maximum and minimum singular values will be closer to 1 in each iteration. Since the maximum value of this singular value is 1, then the contractive matrix will tend to be a unit matrix too. □ So, according to the lemma 4 or the three LMIs Equations (16), (17) and (27), in an iterative loop, the resulted X 0 will be a unitary matrix. By these formulations and constraints, we can conclude the problem of computing the right NCFs in the following theorem; according to this theorem, a brief algorithm will be presented.

Main theorem
An optimal normalized right coprime factors perturbed model Equation (8) This is such that is the maximum difference value between the calculated model and the frequency response set based on the υ-gap metric.
Proof. At each frequency of j , in fact, the problem Equation (10) should be solved such that constraint in Equation (15) is satisfied for unitary matrix X 0 . According to the lemma 1, Equation (13) guarantees Equation (11) and by applying lemma 3, the other two LMI constraints, Equations (16) and (17), are added such that the matrix X 0 will be normalized. According to Lemma 4, the constraint [X * X X * 0 X 0 I ] ≽ 0 is added to the SDP, to guarantee convergence of optimization algorithm. Finally, because of P(s), Q(s) ∈  ∞ , the last two constraints are added.
Except for the condition in [ Δ P Δ Q ] ≤ 1 for uncertainty matrices, there is not any special condition or assumption for matrices and the transfer function matrix. the formulation can be rewritten in the same manner and in Equation (12), instead of X, the variable Y is as follows

Calculation algorithm
There are several software tools like YALMIP [21] for solving the SDP problem resulted in Equation (30). At each frequency, this SDP is solved, and X 0 ( j ) and W( j ) are computed.
In the following, we summarize the proposed method for computing frequency responses of NCFs, in six steps.
1. Initialization: All variables and required settings for the SDP algorithm and also, the frequency index are as j = 1; 2. Determine the initialX at each frequency. As the best initial value, the average value of The algorithm will be finished if the frequency index reached N. • In each frequency, the available samples (K) can be varied.
• The calculation algorithm is independent of the total frequency sample (N) and these samples can be increased, especially in some regions in which the system has a high variation. • After the computation of the frequency responses of the NCFs and weighting matrices, they can be identified by identification techniques.

SIMULATION RESULTS
Demonstrating the ability of this algorithm, in this section, we analyse two MIMO systems as examples.   Example 2. The second example, selected from [18], is a distillation column system with pure delay and ill-conditioned and a nearly singular transfer function matrix. In this example, frequency responses in six operating conditions are used for identification. Figure 4 shows Bode magnitude diagrams of the transfer function matrix of the average frequency response and the identified nominal model by NCFs. Also √ of the model identified in [18] and the computed model by our algorithm are compared together, as shown in Figure 5, which are very close to each other. Of course, in this reference, several optimizations have been done for computing the final optimized model; however, in our algorithm, the optimized NCFs perturbed model is derived directly and by one SDP optimization process.

CONCLUSION
In this paper, the perturbed model of MIMO systems have been computed with frequency responses based on the υ-gap metric minimization. This computation has been performed by introducing a convex optimization problem subjected to some LMI constraints (SDP). Direct computation of the right or left NCFs without additional constraints or conditions can be regarded as the best advantage of this type of formulation; so by having frequency responses, the perturbed coprime factors model of the system and its nominal model can be identified easily. Finally, the extracted algorithm of this formulation has been used for two examples of two references. The results of the perturbed model computed by new algorithm have been compared with those in the references. As shown in the simulated results, our proposed method could directly and optimally compute the frequency responses of NCFs, such that the derived model had the minimum difference with the frequency responses according to the υ-gap metric criterion.
Since, in the presented method and algorithm, any special condition or assumption has not been considered for the transfer function matrix, this method can be used for a vast domain of system. Also, this model can be transformed into another type of perturbed model. The advantage of this coprime decoupling is that Q 1 is diagonal and real and its diagonal elements are non-zero and positive (Q 1 ≻ 0).