A Non-Iterative Approach to Direct Data-Driven Control Design of MIMO LTI Systems

This paper proposes a non-iterative direct data-driven technique that deals with linear time-invariant (LTI) controller design by directly identifying the controller from input-output data without using plant identification. We define the feasible controller parameter set and formulate the problem of designing a controller to match the behaviour of a given reference model as an equivalent set-membership errors-in-variables problem. Then we design the controller parameters by applying recent results in set-membership errors-in-variables identification. Finally, we analyse the effectiveness of the presented technique through simulation examples and experimental results.


I. INTRODUCTION
In recent years, the control community has devoted significant research efforts to the direct data-driven control (DDDC) theory, which is particularly interesting in real-world applications where an accurate model of the plant under control is unavailable.This approach uses input-output data experimentally collected from the plant to design the controller directly.
Active researchers have contributed to several results in this field during the last few years.Campi et al. [5] propose a virtual reference feedback tuning (VRFT) approach.In this one-shot direct data-driven method, a virtual reference signal is the basis of the data-based procedure proposed to design the controller.More precisely, given a set of measured input-output data (u(t), y(t)) and a reference model M (z), one can compute a virtual reference signal r(t) such that r(t) = M (z) −1 y(t).Then, if we feed the closed-loop system and the reference model M (z) with r(t), we can identify a set of controller C(θ, z −1 ) by forcing the closed-loop The associate editor coordinating the review of this manuscript and approving it for publication was Jie Tang .system to match the transfer function M (z).Unlike the standard VRFT for SISO systems, the authors of [19] use variance weighting to achieve a consistent controller estimate with a single input-output data set.Extensions and improvements of this method have been presented in some papers (e.g., [3], [6], [38], [45]).Work [23] presents an alternative data-based approach for controller design called Iterative feedback tuning (IFT).IFT is a data-driven control scheme involving the iterative optimization of a fixed structure controller, tuning its parameters according to an estimated gradient of a control performance criterion.A more recent paper [22] extends the standard IFT approach to the case of multivariable linear time-invariant systems.The Iterative Correlation-based Tuning method (ICbT), proposed in [25], is a data-driven control method in which the controller parameters are tuned iteratively to decorrelate the closed-loop output error, between the designed and achieved closedloop system, from an external reference signal.Although such iterative approaches provide effective techniques to solve the data-driven controller design problem, they show a significant drawback regarding the computational complexity.In fact, several experiments need to be performed at each iteration to solve the gradient estimation problem; furthermore, at least in principle, convergence may require a large number of iterations.In particular, the author of paper [22] shows that for MIMO systems, 1 + n u × n y experiments are required at each iteration, where n u and n y are the number of plant inputs and outputs, respectively.
Direct design of the controller when the collected experimental data are affected by noise is one of the main challenges in the context of DDDC methods.The VRFT method addresses the problem through an extended instrumental variable (IV) approach in work [19].Karimi et al. present an alternative non-iterative approach called Non-iterative Correlation based Tuning (NCbT) in paper [26].This controller tuning approach leads to formulating a controller identification problem where the input is affected by noise while the output is noiseless.Few other solutions have been proposed to cope with the effect of measurement noise (see, e.g., the comparison presented in [43]).
The literature offers several contributions to the complex problem of designing a data-driven controller when the unknown plant shows non-minimum phase (NMP) zeros.Suppose such zeros are not present in the desired reference model.In that case, the closed-loop system is not internally stable because the designed controller leads to unstable zero-pole cancellations (see, e.g., [2] and [40]).
Lecchini and Gevers [32] propose a method, in the context of IFT, to overcome the problem of NMP zero by introducing a flexible reference model.This model has the same poles as the desired reference model, while the parameters of the numerator are free.Inspired by the same idea, [4] uses the flexible reference model in the context of the VRFT approach.The authors of [17] have recently extended this method to multi-input-multi-output (MIMO) plants.
Several other approaches to direct data-driven controller design are available in the literature, although less strictly related to the formulation and solution considered in this work.The interested reader can find a thorough review of most of the available approaches in the survey paper by Wang and Hou [24].
In this study, we present a novel non-iterative approach to direct data-driven controller design, which only relies on a set of input-output experimentally collected data corrupted by bounded noise.The contribution of the work is a substantial extension of our previous conference paper [12], which introduces a set-membership-based direct data-driven controller design technique for stable, minimum-phase (MP) singleinput single-output (SISO) systems.The present contribution significantly improves the work in [12].More specifically, (i) we extend the approach in [12] to deal also with the case of MIMO and NMP systems, and (ii) we provide an original numerical procedure to certify the closed-loop stability of the obtained feedback control system.
First, by only assuming that a bound on the uncertainty affecting the data is available, we formulate the problem of designing a controller in order to match the behavior of an assigned reference model in terms of an equivalent Set-Membership Errors-In-Variables (SM-EIV) identification problem.Then, we design the controller parameters by exploiting previous results by the authors in the field of convex relaxations for errors-in-variables identification (see, e.g., [9]).
The proposed approach has many distinctive features over the other existing procedures in the literature.More precisely, the main novelty of our contribution can be summarized as follows: (i) We only assume to know a bound, not necessarily tight, on the uncertainty corrupting the experimentally collected data.We do not require the knowledge of any statistical information on the noise.Therefore, the results presented in the paper can be successfully applied in those applications where few experimental input-output data are available.Furthermore, we can account for both systematic and random errors, provided that they are bounded: a condition always satisfied in practice.(ii) In contrast to IFT and ICbT approaches that require iterative procedures involving several experiments, the proposed method leads to a one-shot, non-iterative algorithm based on a single experiment.(iii) Unlike the standard VRFT and NCbT approaches, we do not constrain the controller transfer function to depend linearly on the parameters to be tuned.This feature gives the user much flexibility in selecting the desired controller structure.(iv) We propose a unified approach to accommodate both diagonal and non-diagonal MIMO reference models.(v) We propose an original approach to deal with the case of NMP systems.(vi) To the best of our knowledge, we provide the first non-asymptotic data-driven closed-loop stability certification procedure because we only require a finite amount of data.We have organized the paper as follows.In Section II, we formulate the problem, while Section III presents the proposed approach.Section IV presents DDDC for NMP systems, while Section V discusses the stability of the obtained closed-loop system.We show the effectiveness of the presented method in Sections VI and VII through simulation examples and experimental tests, respectively.Concluding remarks end the paper.

II. PROBLEM SETTING A. SISO SYSTEMS
Let us consider the discrete-time linear-time invariant (LTI) single-input single-output (SISO) feedback control scheme depicted in Fig. 1, where q −1 denotes the standard backward shift operator, G(q −1 ) is a stable and MP plant transfer function, K (ρ, q −1 ) is the controller transfer function, ρ is the vector of controller parameters, and M (q −1 ) is the transfer function of a suitable given reference model describing the desired behavior of the controlled plant.
The problem of designing the reference transfer function M (q −1 ) is outside the scope of this work.In paper [14], we propose a procedure to build a reference model guaranteed to fulfil several quantitative performance requirements for both minimum phase (MP) and non-minimum phase (NMP) stable plant by solving a suitably formulated fictitious H ∞ control problem.More precisely, by adequately selecting the weighting functions of the fictitious H ∞ control problem, we can incorporate in the design of M (q −1 ) any given set of specific quantitative performance requirements on robustness, reference tracking, and disturbances rejection.The reader can find a detailed discussion and explanatory simulation examples in [14].The objective of the contribution is to propose an algorithm to design the transfer function of the LTI controller K (ρ, q −1 ) such that the closed loop transfer function T wr (q −1 ) given by T wr (q −1 ) = K (ρ, q −1 )G(q −1 ) 1 + K (ρ, q −1 )G(q −1 ) (1) matches, as close as possible, in some sense, M (q −1 ).This is pursued under the assumption that the plant transfer function G(q −1 ) is unknown, and only a set of input-output data collected by performing suitable experiments on the plant is available.
Let us now introduce the following definitions.

Definition 1 (Model Matching Error Transfer Function):
The model matching error transfer function E(ρ, q −1 ) is defined as the difference between the reference model and the achieved closed-loop transfer function, i.e.

E(ρ, q
Definition 2 (Output Matching Error): The output matching error ϵ(ρ, t) is defined as the signal obtained by multiplying both sides of Equation ( 2) by a reference signal r(t), i.e.
To simplify notation, in the rest of the manuscript, we drop the backward shift operator q −1 from equations and corresponding block diagrams.
Since the output matching error ϵ(t, ρ) in Equation (3) still depends on the unknown plant G, we derive an alternative way of designing the controller K (ρ).For this purpose, we introduce the following result.
Proof of Result 1: Equivalence between (ii) and (iii) in Result 1 follows from the fact that, by imposing ϵ(ρ, t) = 0, the equation can be derived from Equation (3) by means of simple algebraic manipulations.Condition (iii) is now obtained thanks to the fact that GKr(t) = KGr(t) = Kw(t), where w(t) is the plant output sequence obtained by applying the signal r(t) to the plant input (see the block diagram description of the output matching error in Fig. 2).Equivalence between (i) and (ii) is based on the fact that the output of a system is identically zero for all the possible inputs if and only if the system transfer function is identically zero.Remark 1: Result 1 plays a crucial role since it suggests a way for turning condition (i) on the model matching error, which depends on the unknown plant transfer function G, into condition (iii) which, on the contrary, depends only on the output sequence w(t) collected by applying the signal r(t) to the plant (see Equation ( 6)).
Remark 2: Condition (ii) considered in Result 1 can be approximated in practice by the condition ϵ(ρ, t) = 0 for a reference signal r(t), which is persistently exciting in the sense that its frequency spectrum is rich enough to properly excite the dynamics of both M and G.For more details see [12].
Remark 3: It is worth noting that several previous papers on the subject (see, e.g., [19], [42], and [46]) introduce, at this stage, the following approximation: The meaning of such an approximation is that the actual sensitivity of the controlled system to be obtained with the designed controller will be equal to the ideal sensitivity (1 − M ), a condition which, in turn, implies perfect matching of the input-output reference model.Such an assumption, which was also made in our preliminary conference paper ( [12]), is no more required in this work.In fact, by exploiting the set-membership technique presented in the next sections, the controller K (ρ) in Equation ( 6) can be directly designed from the I/O collected data without assuming that the final controlled system transfer function will perfectly match the desired reference model (a condition which cannot be guaranteed in the general case).

B. MIMO SYSTEMS
In this section, we consider the discrete-time linear-time invariant (LTI) feedback control scheme shown in Fig. 1 as a multi-input/multi-output (MIMO) system where G ∈ C n×n , K(ρ) ∈ C n×n and M ∈ C n×n are linear discrete-time transfer function matrices of a stable and MP plant, a controller and a given reference model respectively.As far as the selection of the reference model is concerned, the same considerations discussed in Section II-A for SISO systems apply.In this work, we assume that the plant to be controlled is a square system in the sense that the input and the output vectors, u(t) and w(t) respectively, have the same dimension n.We make this assumption in order to simplify the presentation for the sake of readability.However, all the presented results and algorithms also apply to non-square systems.In this paper G ij , M ij and k ij denotes the (i, j)-element of G, M and K matrices, respectively.We introduce the following result.
Result 2: The following three conditions are equivalent where, I is the identity matrix of appropriate dimension (i.e., n × n).Proof of Result 2: Equivalence between (i) and (ii) is based on the fact that the output of a system is identically zero for all the possible inputs if and only if the system transfer function is identically zero.
The equivalence between (ii) and (iii) in Result 2 can be proved by noting that condition (ii) is equivalent to Let us now introduce a generic signal χ(t) ∈ R n such that Since Equation ( 12) holds for all r(t) ∈ R n , then it holds, in particular, for the signal χ(t) ∈ R n satisfying (13).As a consequence, we get

III. SET-MEMBERSHIP-BASED DIRECT DATA-DRIVEN CONTROL TUNING (SM-DDDC)
In this section, we propose a set-membership-based algorithm for DDDC design for both SISO and MIMO systems.

A. SM-DDDC FOR SISO SYSTEMS
Here we assume that a set of N input-output plant data are collected experimentally by applying a suitable (persistently exciting) signal r(t) to the plant.The output measurements y(t) are assumed to be corrupted by bounded additive noise according to where is the noiseless output of the system, while the noise η(t) is assumed to range within a given bound η as follows Substitution of Equation ( 15) into ( 6) leads to the following equation Because of Result 1, equation ( 18) represents the condition to be satisfied by the controller K (ρ) to solve the model matching problem described in Section II-A.Since such a condition depends on the uncertain variable η(t), equation ( 18) does not have an exact unique solution.Therefore, according to the set-membership paradigm (see, e.g., [9]), we consider the set of all LTI controllers K (ρ) of given order n ρ which are the solutions to equation (18) corresponding to all the possible noise sequences bounded as in (17).We call this set the Feasible Controller Set (FCS).Practically speaking, the FCS represents the solution space of our problem from which we have to pick the specific controller to be implemented on the plant.Fundamental properties of FCS are discussed in Section III-C, while in Section III-D, we propose a way to select optimally from the FCS the controller to be implemented.The FCS for SISO systems is formally defined as follows.
Definition 3 (Feasible Controller Set for SISO Systems): The Feasible controller set for SISO systems is defined as the set of all the controllers belonging to the following given class such that the equation is fulfilled for at least one noise sequence η(t) satisfying the bound |η(t)| ≤ η, ∀t = 1, . . ., N .
121674 VOLUME 11, 2023 Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.
The Feasible controller set can then be described as It is worth noting that the controller class K (n a ,n b ) is general in the context of linear control since it includes all the linear time-invariant (LTI) controller of order at most max(n a , n b ).
Since the controller class is parametrized by ρ, we can conveniently replace the set D K with the Feasible controller parameter set defined as follows.
Definition 4 (Feasible Controller Parameter Set for SISO Systems): The feasible controller parameter set is the set of all the controller parameters ρ such that K (ρ) ∈ D K .
Since Equation ( 21) not only depends on ρ but also on the uncertain variable η(t), we consider the following set.
Definition 5 (Extended Feasible Controller Parameter Set for SISO Systems): The extended feasible controller parameter set D ρ is the set of all the controller parameters ρ and noise sequences η(t), t = 1, . . ., N , such that Based on the definition of the extended feasible controller parameter set, we can formulate the following result.

Result 3 (Structure of the Extended Feasible Controller Parameter Set for SISO Systems):
The extended feasible controller parameter set can be written in the following form where Proof of Result 3: From definition 5, D ρ is the set of parameter ρ ∈ R n ρ such that the following condition holds which, in turn, is equivalent to where L(t) is defined as in (23).The statement of Result 3 finally follows from the fact that

B. SM-DDDC FOR MIMO SYSTEMS
In this Section, a set-membership algorithm for DDDC design is proposed to deal with the case of MIMO LTI plants.We assume that a set of N input-output plant data are collected through an open-loop experiment by applying a suitable (i.e., persistently exciting) reference signal r(t) to the plant G.The output measurement y(t) is assumed to be corrupted by bounded additive noise according to where the noise variables η(t) acting on the noiseless output w(t) are assumed to range within given bound η , that is Substitution of Equation ( 27) into (11c), leads to the following equation where, due to the presence of noise on the output measurements, the controller design depends on the uncertain variables η(t) as in the SISO case.Now, let us introduce the Feasible Controller Set for MIMO systems defined as follows.
Definition 6 (Feasible Controller Set for MIMO System): The feasible controller set for MIMO system is the set of all the controllers belonging to the following given class where the generic entry of such a matrix is described as, such that the equation is fulfilled for a noise sequence η(t) satisfying the bound The Feasible controller set can then be described as Since the controller class considered here is parametrized by ρ, we can conveniently replace the set D K with the Feasible controller parameter set defined as follows, Definition 7 (Feasible Controller Parameter Set for MIMO System): The feasible controller parameter set D ρ is the set of all the controller parameter vectors ρ such that K(ρ) ∈ D K .
Since Equation ( 21) not only depends on ρ but also on the uncertain variables η(t) and χ(t), we consider the following set.
Definition 8 (Extended Feasible Controller Parameter Set for MIMO Systems): The extended feasible controller parameter set D ρ is the set of all the controller parameter vectors ρ, noise sequences η(t), t = 1 . . ., N and uncertain signal samples χ(t), t = 1, . . ., N such that Based on the definition of the Extended Feasible controller parameter set, we can formulate the following result.
Result 4 (Structure of the Extended Feasible Controller Parameter Set for MIMO System): The extended feasible controller parameter set can be written in the equivalent form by introducing some additional variables, Z ij and Q ij , called partial outputs (see Appendix A and [11] for more details) Proof of Result 4: See Appendix A. Remark 4: In many real-world applications, the offdiagonal transfer functions of the controller matrix K are tuned in such a way that the interactions between the outputs are decoupled, and therefore G is assumed to be diagonalizable by output feedback (see, e.g., [34]).Thanks to Result 4, this decoupling can be satisfied by choosing a diagonal reference model M. In this case, the formulation of the extended feasible controller parameter set enjoys additional properties that can be exploited to reduce the number of slack variables involved in its description, and thus reduce the computational effort required to compute the controller.See Appendix B for a detailed mathematical description.For more details about controller decoupling in linear multivariable systems, see, e.g., [34], [39], and [44].

C. FUNDAMENTAL PROPERTIES OF THE FEASIBLE CONTROLLER SET
The importance of the feasible controller set comes from the fact that D K is the set of all the controllers having order less than or equal to max(n a , n b ), which are consistent with the assumption that the output matching error is identically zero for at least one of the possible feasible noise sequences.In turn, this implies that no controller in the considered class K (n a ,n b ) can solve the problem if and only if the feasible controller set D K is empty.If that is the case, we must enrich the model class by considering higher values for n a and/or n b (i.e., higher controller order).The converse is not true in general, i.e., if D K is not empty, it is not guaranteed to contain the ideal controller defined by which exactly solves the model-matching problem.To ensure this desirable property, we also need to guarantee that the controller class is rich enough.Formally we introduce the following result.

Result 5 (Emptiness of the Feasible Controller Set
Then, for all N ≥ 0, the feasible controller set is not empty and always contains the ideal controller K id .
Proof of Result 5: The above result immediately follows from the definition of the ideal controller.
If n a , n b < n G + n M , either the ideal controller performs cancellations with the plant, and the feasible controller set is non-empty for all N ≥ 0, or the feasible controller set does not contain the ideal controller and eventually could become empty when N increases.
Remark 5: Since we use physical sensors to collect the output data y experimentally, the measurements are corrupted by the effect of disturbances and errors.Our approach models the uncertainty affecting the collected data with the unknown but bounded variable η.Since η enters into the definition of the feasible controller set, the FCS size could become significant in the presence of large uncertainty η.In such a situation, assuming that a sensor model is available, we can design an interval observer around the sensor to obtain an improved estimate of the measured variables, thus reducing the size of the uncertainty η.The design of such an observer is outside the scope of the paper.The interested reader is referred to the survey papers [27], [28].

D. CENTRAL ESTIMATE
Given the set D ρ of all the feasible controller parameters ρ, we pick a specific single value in the set to design the controller K to be actually implemented in the feedback control system of Fig. 1.According to the set-membership identification theory, we select a single point in D ρ by looking for the value of the parameter ρ that minimizes the distance in the ℓ ∞ norm from the farthest point in the feasible controller parameter set, i.e.
The estimate ρ c computed by solving (35)  Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.
the center of the minimum-volume-box outer-bounding D ρ and can be computed by exploiting the convex relaxation approach proposed by [9].In particular, for each single component ρ x of the parameter vector ρ the central estimate ρ c x is given by where Controller parameter bounds ρ x and ρ x implicitly define the Controller parameter uncertainty intervals (CPUIs) given by Problem ( 37) is, in general, a hard non-convex optimization problem that falls into the class of semi-algebraic optimization problems, widely studied in recent years.More specifically, it has been shown that the global optimum of a constrained semi-algebraic program can be approximated arbitrarily well by exploiting either the moment-basedapproach [29] or the sum-of-squares-based decomposition approach proposed in [15] and [35].Methodologies proposed in [15], [29], and [35] allow the user to construct a sequence of convex linear matrix inequality (LMI) problems, guaranteed to converge to the global optimum of the original non-convex polynomial problem as the order of relaxation goes to infinity (see the book [30] and the references therein for details).However, the direct application of such methods to large-scale identification problems (i.e., a large number of parameters to be estimated and/or a large set of experimental input-output data) might lead to intractable LMI problems due to the required memory storage and/or computational time.To overcome this limitation, ad hoc approaches have been proposed in [7], [8], [9], and [10], aimed to reduce the computational complexity by exploiting some structural features of the polynomial optimization problems arising from the context of system identification.In particular, it is possible to show that problem (37) enjoys the same sparsity structure as the problem considered in [7], [8], [9], and [10].Thus computationally effective implementation can be applied to solve DDDC problems with several hundreds of input-output data.
Remark 6: We remark that in the proposed work, we need to solve a non-convex optimization problem.This fact is a major drawback compared with the most common DDDC methods, such as VRFT and NCbT, where the controller design only requires the solution of convex optimization problems.However, in these methods, the user is asked to design a linear controller such that the denominator of K (ρ) is fixed (usually 1 − q −1 ), which has a big limitation on the controller structure that needs to be designed (see, e.g., [38] and simulation Example 1).Moreover, this limitation can affect the output results for NMP systems, as discussed in the next section and illustrated in the simulation Example 3.
To overcome this problem, few methods in the literature have been used to extend the basic SISO VRFT and NCbT to obtain a fully parameterized controller (see, e.g., [38]).However, these methods failed to have good results when the data is corrupted by noise and to design a fully parameterized controller for MIMO systems.Therefore, one of the main distinctive features of the SM-DDDC approach is that the controller transfer function (both for SISO and MIMO systems) does not need to depend linearly on the parameters to be tuned, and no statistical information about the noise is assumed to be a-priori available.

IV. DIRECT DATA-DRIVEN DESIGN FOR MIMO NON-MINIMUM PHASE PLANT
An important limitation related to the internal closed-loop stability arises when the plant possesses transmission zeros outside the unit circle.Pole-zero cancellation issues may occur when applying model-matching design techniques to NMP systems.This section studies a data-driven design methodology to handle the case of NMP systems.
As stated in the SISO case in [13], when the plant has NMP zeros that are not included in the reference model, the model matching controller may lead to instability.In fact, in this case, the ideal controller K id (q −1 ) will certainly include in its denominator all the transmission zeros of G that are not included in the reference model M.However, unstable poles may arise in the ideal model matching controller and also in the case of minimum phase plant, as shown in [13].Moreover, the notion of NMP zero is much more complex for MIMO systems than for the SISO case.In particular, to discuss the internal stability of NMP MIMO systems we need to refer to the notion of transmission zeros and their direction.On this basis, and for the self-consistency of the paper, the following definitions for zeros for the multivariable system are introduced.
Definition 9 (Transmission Zeros [31]): The finite transmission zeros polynomial z(q) of a system G(q) is the greatest common divisor of all numerators of all minors of order n r of G(q), where n r is the normal rank of G(q).Definition 10 (Transmission Zero Direction [31]): If G(q) possesses a transmission zero at q = z i , then there exist non-null vectors known as the zero input direction u z i and zero output direction y z i .These vectors satisfy the following conditions: u H z i u z i = 1, y H z i y z i = 1, and furthermore, they fulfill the relationships G(z i )u z i = 0 and y H z i G(z i ) = 0.In principle, we may obtain u z and y z from an SVD of G(z) = U V H ; and we have that u z is the last column in V and y z is the last column of U , corresponding to the zero singular value of G(ζ ).For more details, see e.g., [31] and [33].
Result 6 (Internal Stability of MIMO Feedback System [21]): With reference to the block diagram depicted in Fig. 1, if G is stable and has NMP transmission zero at z i with output direction y z i , then the feedback system with controller K will be internally stable if the following interpolation constraint is satisfied: In words, Result 6 states that M must have the same transmission zeros of G in the same output direction.It is important to notice that constraint (39) is a function of the transmission zeros and has no direct relationship with the zeros of the elements of G.
Remark 7: If the reference model M is diagonal (M ij = 0 ∀i ̸ = j), the NMP transmission zeros will be present in each element of M. Thus, for a diagonal structure of the reference model, one does not need to be concerned with the zero output direction since the constraint in the Result 6 will be satisfied because T wr (z i ) = 0.For more details, see, e.g., [16].
To detect and locate the transmission zeros and their output direction, we introduce the notion of prospective plant for MIMO systems.
Definition 11 (Prospective plant G * (ρ)): Given a reference model with transfer function M and a controller with transfer functions K(ρ), the prospective plant is defined as the plant estimate given by the following transfer function Based on the notion of the prospective plant, we propose the following Two-Stage design in Algorithm (1), which allows the user to detect if the designed controller could lead to unstable pole-zero cancellation.

Non-diagonal reference model
2.2.1.Add all the NMP transmission zeros of K (ρ) and their directions to the numerator of the M; 2.2.2.Modify the NMP transmission zero of M to have an output direction equal to the process and with a different input direction; 2.2.3.Move the effect of the NMP zero to a specific output using e.g., block triangular structure or another suitable method (for more details see [16]); 2.2.4.Compute K (ρ), using the modified M, as in Section III.
Remark 8: It is worth noting that, in case the prospective plant G * (ρ) shows the unstable poles of the controller as zeros, this suggests that the unstable poles appeared in the controller to match the selected reference model M. Therefore, the chosen reference model has to be modified.On the contrary, if G * (ρ) doesn't show unstable zeros corresponding to the unstable poles in the designed controller K, this suggests that the selected reference model M did not induce unstable zero-pole cancellations.Therefore, we don't need to modify the chosen reference model.
Remark 9: The main distinctive features of the proposed approach (two-stage design) with respect to those already available in the literature for NMP systems (e.g., [1], [4], [17], [36], and [41]) are as follows: (i) no a-priori information on the NMP zeros location is needed; (ii) the proposed strategy does not involve the optimization of the reference model.

V. DATA-DRIVEN STABILITY CERTIFICATION
In this section, we propose a data-driven approach, in the same spirit as the controller certification paradigm discussed in [20], to check if the designed controller provides closedloop stability.In particular, according to the approach proposed in [42], we assume that a stabilizing controller K s is a-priori known.Such a controller only provides stability of the closed-loop system; it is not expected to solve the model-matching problem formulated in Section II.Then, we consider the matrix transfer function defined by where M s .= (I + K s G) −1 GK s is the complementary sensitivity function obtained with the given stabilizing controller K s .The following result provides a sufficient condition for closed-loop stability.
Result 7: [Closed-Loop Stability] The closed-loop system obtained with the designed central controller K(ρ c ) is stable if (ρ c , z) is BIBO stable and where (ρ c , z) denotes the transfer function of the system (ρ c , q −1 ) in the Z-domain.∥ • ∥ ∞ and σ x (•) denote the H ∞ norm of a system and the maximum singular value of a matrix, respectively.This result directly follows from the application of the small-gain theorem and it is an extension to the MIMO case of the result presented in the work [42], where the reader can find a detailed discussion and proof.Moreover, as shown in [42], stability of (ρ c , z) is ensured any time K(ρ c ) is BIBO or contains the same number of poles in z = 1 as (I − M s ) −1 .
Result 7 cannot be applied in practice in its original form (42) since (ρ c , e iω ) depends on G which is unknown.In order to overcome this problem, we propose here a procedure to compute a data-driven upper bound on σ x (ρ c , e iω ) for any ω ∈ , where is a finite subset of [0, ∞).We also discuss how to construct to ensure stability despite the fact that is a finite subset of [0, ∞).
121678 VOLUME 11, 2023 Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.
First, let us describe the (i, j)-entry of the matrix (ρ c , q −1 ) as a discrete-time transfer function of the following form where d is the vector of unknown transfer function parameters.Although the transfer function order n is not exactly known, the following upper bound can be easily derived from the definition of (ρ c , e iω ) given in (41): where n M are the orders of K ij (ρ c ) and [M s ] ij , and n G [i,j] is an upper bound of the true plant dynamic order.We can formally state the data-driven estimation problem to be solved as follows.
Problem 1 (Data-Driven Set-Membership Estimation of σ x (ρ c , e iω ) ): Let (q −1 ) be an unknown matrix transfer function satisfying Equation (41).Assume that: (i) A set of input-output data are experimentally collected on the plant G. (ii) The collected output is corrupted by bounded noise according to equations ( 15)-( 17).(iii) The entries of matrix (q −1 ) are described by the parametric model in Equation (44).Compute an upper bound on σ x (ρ c , e iω ) , ∀ω ∈ .The solution of Problem 1 is provided by the following result.
Result 8: For each fixed value of ω, the solution to Problem 1 can be obtained by computing the optimum of the following optimization problem: (d, e iω ) = Udiag(σ )V H (45b) where A H denotes the conjugate transpose of A, σ is the vector of singular values of (d, e iω ) listed in decreasing order, i.e., σ i+1 < σ i , diag(σ ) is a diagonal matrix with the singular values on the diagonal, and U, V are complex matrices of suitable dimension.The proof of Result 8 is obtained by simply analyzing the meaning of the constraints added to the problem.Constraints (45b) and (45c) ensure that the minimizer (U * , V * , σ * ) defines the singular value decomposition (SVD) of (d * , e iω ).Constraints (45d) and (45e) provide a data-based equivalent of condition (41) obtained by applying the same procedure employed in the reformulation of the extended feasible controller parameter set, as presented in Appendix A.
Problem ( 45) is a polynomial optimization problem since both the cost function and all the constraints are multivariate polynomial functions of decision variables.In particular, constraint (45b) is equivalent to 2n 2 scalar polynomial equality constraints of degree 3. Notice that matrix optimization variables U and V are complex; thus, we need 4n 2 scalar real variables to represent them as real polynomials.
Consequently, Problem 1 can be solved to the global optimum by means of suitable convex relaxation techniques by selecting a sufficiently large order or relaxation; for low values of the order of relaxation, a conservative upper bound is obtained (see the discussion in Section III-D).Applying this procedure for all ω ∈ allows us to over-bound σ x ( (ρ c , e iω )) in the set of frequencies .
The condition σ x ( (ρ c , e iω )) < 1 for all ω ∈ is a reliable approximation of the stability condition in (7) when is sufficiently rich (i.e., the points are selected with a fine grid and in such a way to covers a wide range of frequencies).However, the condition becomes rigorously sufficient only when = [0, ∞), but this leads (both in the presented setting and in the approach considered in [42]) to the formulation of semi-infinite optimization problems (i.e., with an infinite amount of constraints), which are intractable in general.
To overcome this difficulty, we introduce the following assumption: Assumption 1 (Lipschitz Continuity of the Maximum Singular Value): The maximum singular value σ x ( (ρ c , e jω )) is Lipschitz continuous in ω, and an upper bound L b on the Lipschitz constant is available.Under this assumption, ∥ (ρ c , z)∥ ∞ is over-bounded by where δ ω represents the maximum distance between two consecutive frequencies in .Therefore if then K(ρ c ) stabilizes the unknown plant.Remark 10: In general, the constant L b is not available a-priori.However, we remark that this value can always be chosen as a very large conservative value.In this case, we must select the frequency gridding parameter δ ω small enough to maintain the product δ ω L b small enough.
Remark 11: If condition (47) is not satisfied, a possible countermeasure is to redesign the controller using a larger amount of experimentally collected data.Increasing the number of input-output data will reduce the size of the feasible controller parameter set.Consequently, the distance between the central estimate and the ideal controller becomes smaller and, in turn, ∥ (ρ c )∥ ∞ decreases.
Remark 12: Differently from the approach proposed in [42]: Comparison of frequency responses: designed feedback control system with SM-EIV (solid thick-line), NCbT method(solid thin-line), VRFT method(dashed) and reference model (dotted).
From the comparison, we see that the controlled system obtained with the proposed method performs better than the NCbT and VRFT methods both in the frequency and in the time domain.

B. EXAMPLE 2
In this example, we present a comparison between the SM-EIV approach proposed in this paper, the VRFT method proposed in [19] and the method proposed in [18].The plant transfer function used for generating the data is given by while the assigned non-symmetric reference model M 2 is considered for the same plant, given by The input signal for the SM-EIV method is defined as r(t) = [s 1 (t), s 2 (t)] ⊤ , where s i (t) is a random signal uniformly distributed in the range [−1, +1] with length N = 100, for i = 1, 2. The input signal for the VRFT method is u(t) = [A 1 (t), A 2 (t)], where A i (t) is a random signal with length N = 5000, for i = 1, 2. Finally, the input signals for the method proposed in [18] are two random input sequences u 1 and u 2 with length N = 5000 such that, the sequence u 1 is used to feed the first channel, then the same input is switched to the second channel and u 2 is used to separately feed the two channels, analogously to what has been done for u 1 .Note that, for the SM-EIV and the VRFT we need one experiment while for the method proposed in [18] we need n × n experiments, where in this example n = 2.
To establish a fair comparison with the method of [19] and [18] two different tests have been done, one with bounded noise uniformly distributed in the range [− η, + η] and another one with stochastic zero-mean white noise, such that both systems are characterized by SNR w ∼ = 25 dB.Nonetheless, it should be noticed that the quality of the step response of all the previous methods is the same for both tests (stochastic and bounded noise) in terms of overshoot, performance in channel decoupling, settling time and overall shape.
The controller order for the SM-EIV method is selected by trial starting from n = 1 and by increasing n until the feasible set is not empty.Since the feasible set is not empty for n = 1, we select the following controller structure b [11]   0 + b [11]  1 q −1 1 + a [11]  1 q −1 b [12]   0 + b [12]  1 q −1 1 + a [12]  1 q −1 b [21]   0 + b [21]  1 q −1 1 + a [21]  1 q −1 b [22]   0 + b [22]  1 q −1 1 + a [22]  1 While, the controller structure for the VRFT and for the approach proposed in [18] is chosen to be linearly parameterized with a fixed known denominator, as follows b [11]   0 + b [11]  1 q −1 1 − q −1 b [12]   0 + b [12]  1 q −1 1 − q −1 b [21]   0 + b [21]  1 q −1 1 − q −1 b [22]   0 + b [22]  1 q A closed-loop noiseless experiment with the controller given by the proposed approach and the one returned by MIMO VRFT and the one proposed in [18] is illustrated in Figure 5.As can be seen from Fig. 5, the proposed approach (SM-EIV) provides a perfect decoupling, ensuring that the closed-loop system is diagonalized, thanks to the diagonal reference model.The implementation relies on the simplified description of the extended feasible controller parameter set provided in Appendix B. On the contrary, the VRFT method for MIMO systems and the method proposed by [18] don't guarantee decoupling because using a finite amount of data introduces a bias in the controller estimate.

FIGURE 5.
Step responses: designed feedback control system with the SM-EIV approach (solid thick-line), the VRFT method (dashed), the method proposed by [18] (solid thin-line), reference model (dashed-dotted) and reference signals (dotted).Notice that dashed-dotted and solid lines are perfectly overlapped.

C. EXAMPLE 3
In this example, the proposed approach is employed to tune an NMP SISO controller to be compared with the standard NCbT introduced in [42] and the VRFT approach (see, e.g., [19]) when suitable convex stability constraint is imposed (see, e.g., [42]).The plant considered in this example is taken from [37] and has the following transfer function Step responses: designed feedback control system with the SM-EIV approach (solid-thick), the NCbT method (dashed-dotted), the VRFT method (dashed), reference model (dotted) and the modified reference model obtained (solid-thin).Notice that, dashed-dotted and solid-thick lines are perfectly overlapped.
while the assigned reference model has a transfer function The system is excited by a random input signal r(t) uniformly distributed in the range [−1, +1].The plant output signal w(t) is corrupted by a random additive noise η(t), uniformly distributed in the range [− η, + η].The chosen error bound η is such that the signal-to-noise ratio is 26 dB.Input-output samples are collected with a sampling time T s = 0.1 s and, discrete-time models of G(s) and M (s) are obtained, for simulation purposes, through ZOH discretization method, according to [37].
In this example the controller K (ρ) is obtained by computing the central estimate ρ c j = (ρ c j + ρ c j )/2 through the convex relaxation approach proposed in [9] for a relaxation order δ = 2.The general LTI controller structure in (19) is considered here where n a = n b = n and the controller order n = 2 is selected by trial and error starting from n = 1 and by increasing n until the feasible set is not empty.The final selected controller structure is as follows where, by taking the central estimate of the controller parameters the following controller transfer function is obtained ) .
Since the obtained controller has an unstable pole, according to the Two-Stage procedure presented in Section IV, we have to compute the prospective plant G * in Equation ( 40) the equation can be derived, as shown at the bottom of next page.
Since G * has one NMP zero, we can conclude that the unknown plant G is an NMP system and the NMP zero of G * (i.e., the unstable pole in K ) is expected to be a good estimate of the NMP zero of G. Therefore the unstable pole of K (ρ) is added to the reference model M as an NMP zero.Finally, the problem is solved by exploiting the modified reference model with a controller order n = 3.The final obtained controller transfer function (corresponding to the central estimate of the controller parameters) is the following As for VRFT NCbT, a PRBS signal 255 samples with unity amplitude is used as input to the system.Four periods of this signal are used to design the controller, N = 1020.The output is disturbed by a zero-mean white noise such that the signal-to-noise ratio is about 26 dB.The length of the instrumental variable l and the rectangular window l 2 (used for computing the H ∞ norm estimate for the stability constraint) are found to be 10 and 120 respectively, by trial-and-errors.
Fig. 6 displays the comparison between the reference model and the obtained closed-loop system in terms of step responses respectively for SM-EIV, VRFT and NCbT methods.From the comparison, we see that the controlled system obtained with the proposed method performs better than the NCbT and VRFT methods.

D. STABILITY CERTIFICATION EXAMPLE
In this example, we apply the data-driven controller certification method proposed in Section V to the controller obtained in Example 1 (Section VI-A).The central estimate controller to be certified is To apply the proposed procedure, we assume to know the complementary sensitivity function M s achieved by a stabilizing controller K s .As discussed in [42], if the plant is stable and minimum phase, then for any stable M s , there exists a stabilizing ideal controller K s .Moreover, (ρ c , z) is stable thanks to the fact that (1−M s (z)) −1 and K (ρ c , z) both one and have only one pole in z = 1.As a direct consequence of those observations, we select M , defined in Example 1, as a good candidate for M s .Furthermore, we assume the knowledge of an upper bound on the Lipshitz constant of σ x ( (e iω )).As discussed in Remark 10, estimating this quantity is hard.However, any sufficiently large quantity can be safely assumed if the upper bound on σ x ( (e iω )) is computed for a sufficiently large set of frequencies.Accordingly, we select 15 logarithmically equispaced frequencies between 1 × 10 −2 rad s −1 and 2.154 rad s −1 , and we assume that the Lipshitz constant of σ x ( (e iω )) is at most 200 dB/dec.
We evaluate the solution to (45) for all selected frequencies using the same data used in Example 1 to estimate the 121682 VOLUME 11, 2023 Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.controller parameters.Then, using the information about the Lipshitz constant, we get a bound for all frequencies not considered in the optimization.Figure 7 shows the obtained results.According to Equation ( 46), we get that the upper bound on σ x ( (e iω )) is −9.8 dB ≈ 0.3236.Since this value is well below 1, we conclude that K c is a stabilizing controller for the unknown plant.
We remark that the computed upper bounds on σ x ( (e iω )) can be made less conservative by increasing the number of samples used in the formulation of optimization problem (45), at the cost of increased computational complexity.Moreover, checking the stability of controllers provided in Example 2 and Example 3 consists of the same procedure, provided that some appropriate M s and L b are selected.

VII. EXPERIMENTAL RESULTS AND DISCUSSION
The algorithm presented in Section III, has also been tested on the experimental input-output data collected on a test bench MIMO electronic filter with 2 inputs and 2 outputs.Fig. 8 shows the experimental setup used to collect the data.
The system structure is reported in the block diagram depicted in Fig. 9, where G 11 is the transfer function of a second-order low-pass filter with two complex-conjugate poles, characterized by a natural frequency of 95 Hz and a damping factor of 0.6.The transfer function has been practically built in the form of a Sallen-Key circuit.G 21 is a transfer function of a third-order low-pass filter with a couple of complex conjugated poles with a natural frequency 120 Hz and a damping factor of 0.5, and a real pole at 160 Hz.The physical realization was done by means of a Sallen-Key circuit implementing the complex conjugated poles pair, and an RC circuit, for the additional real pole.G 12 is the transfer function of a high-pass filter with a pair of complex conjugated zeros with a natural frequency of 17 Hz and damping factor of 0.2, and a pair of complex conjugated poles with a natural frequency of 83Hz and a damping factor of 0.5.The transfer function was implemented by means of a Tow-Thomas circuit.G 22 is the transfer function of a first order low-pass filter with a real pole at 80 Hz built in the form of a standard RC filter.We point out that, in this example, we do not assume any a-priori information on the plant G. G * (q −1 ) = 0.087901(1 + 0.964q −1 )(1 − 1.0515q −1 )(1 − 0.2495q −1 ) (1 + 0.9527q −1 )(1 − 0.9189q −1 )(1 − 0.8912q −1 )(1 − 0.2384q −1 ) In order to achieve reasonable tracking and decoupling between y 1 (t) and y 2 (t), we select the following discrete-time reference model The system has been excited with the input signal, r(t) = [s 1 (t), s 2 (t)] ⊤ , where s i (t) is a random sequence of 200 samples, uniformly distributed within the range [−1, +1]V, for i = 1, 2. A National Instruments PXI equipped with a NI-6221 DAQ board has been used to generate the input signals s i (t) and to collect the signals r(t) and y(t) at a sample rate of 4 kHz.The upper bound on the measurement errors is derived from the precision of the measurement equipment which is given by η = 8 mV.The software sparsePOP and MOSEK has been used to solve the underlying optimization problems.The central estimate of the controller parameters for the SM-EIV method (ρ c SM −EIV ) has been obtained by exploiting the convex relaxation approach proposed in [9] for a relaxation order δ = 2.
The controller order for the SM-EIV method is selected by trial starting from n = 1 and by increasing n until the feasible set is not empty.Since the feasible set is not empty for n = 1, we select the following structure for the controller b [11]   0 + b [11]  1 q −1 1 + a [11]  1 q −1 b [12]   0 + b [12]  1 q −1 1 + a [12]  1 q −1 b [21]   0 + b [21]  1 q −1 1 + a [21]  1 q −1 b [22]   0 + b [22]  1 q −1 1 + a [22]  1 q −1        The final controller parameters for the SM-EIV (ρ c SM −EIV ) method are reported in Table 1.As far as the parameters value of a [ij] 1 , ∀i, j = 1, 2 are concerned in the SM-EIV technique, the computed central estimated for these parameters are equal to −1.A comparison between the reference model and the obtained closed-loop system for the SM-EIV method setting the value of the parameter to the central estimate ρ c SM −EIV ) in terms of the step response is presented in Fig. 10, from which we see, using 200 samples data corrupted by noise, that the output of the controlled response for the SM-EIV approach and the desired output from the reference model are practically indistinguishable.

FIGURE 10.
Step responses: feedback control system with the SM-EIV approach (solid), reference model (dashed) and reference signals (dashed-dotted).solid lines are perfectly overlapped.

VIII. CONCLUSION
In this work, we have proposed an original approach to the problem of designing linear multiple-input multipleoutput (MIMO) controllers directly from a set of input-output data experimentally collected on the plant to be controlled.Assuming that a bounded noise affects the output measurements, we have formulated the controller design problem as a peculiar input-error set-membership identification problem solved by adapting previous results on errors-in-variables identification proposed by some of the authors of this work.The presented method has many advantages over previously proposed non-iterative direct data-driven design algorithms.Firstly, we are not bounded to consider linearly parametrized controller structures.Secondly, we can deal with diagonal and non-diagonal MIMO reference models, and thirdly, we do not need prior information on the location of possible non-minimum phase transmission zeroes.The proposed approach can be extended straightforwardly to non-square systems.Thanks to these distinctive features, the approach considered here is significantly more flexible than the previously available ones.We have shown the effectiveness of the proposed approach through simulation examples and the application to a laboratory test bench.
121684 VOLUME 11, 2023 Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.

FIGURE 1 .
FIGURE 1.Feedback control system to be designed compared with the reference model M(q −1 ).

Result 1 :
The following three conditions are equivalent

FIGURE 2 .
FIGURE 2. A block diagram description of the output matching error ϵ(t , ρ).

FIGURE 7 .
FIGURE 7. Bode plot of σ x (ω): true value computed using the unknown plant (solid line), data-driven estimated upper bound (dashed), and point evaluated though the solution of optimization problems (45) (asterisks).

FIGURE 8 .
FIGURE 8.The experimental MIMO system used as a test bench.

FIGURE 9 .
FIGURE 9. Block-diagram description of the MIMO circuit considered in the experimental test bench section.