On PID Design Constraints in Relation to Control of Strictly Metzler Linear MIMO Systems

Abstract

For time-invariant Metzler linear MIMO systems, this paper proposes an original approach reflecting necessary matching conditions, specifically structural system constraints and necessary positiveness in solving the problem of MIMO PID control. Covering the matching conditions by the supporting structure of measurement, refining the controller and system parameter constraints and introducing enhanced equivalent system descriptions, the reformulated design task is consistent with PID control law parameter representation and is formulated as a linear matrix inequality feasibility problem. Characterization of the PID control law parameters is permitted to highlight dynamical properties of the closed-loop system and the structural influence of the control derivative gain value in the design step. For the first time, the paper comprehensively sets the synthesis standard for PID control of MIMO Metzler systems because others for the given task have not been created at present.

1. Introduction

Representing the system state and variable positiveness of systems in different domains [1], positive systems act as a specific class of systems of technical processes. To cover the class of plants using the model with non-negative parameters, the distinctive connections with the Metzler structure of system matrices [2] for notational simplicity implies a rather common notation of them as Metzler systems. The main areas of applicability are switched systems [3] and multi agent systems [4]. Since the existence of controllers that stabilize this class of plants is a distinctive problem, suitable publications are focused on the sophisticated techniques applicable for representation of different positive constraints [5,6]. To achieve the necessary closed-loop system state positivity with respect to parameter boundaries and then semi-definite programming [7], the implementation of non-symmetrical bounds [8] as well as the combined linear programming method [9] have been proposed, but the specific problems concerned with solver interactions and parameter constraints remained open.

Tractable ways of including a PID controller into control structures have appeared, where stabilizing tasks formulated with the inclusion of additional performance requirements and constraints are also of interest [10]. The resulting closed-loop system provides stable system responses and, if the design is covered by adequate matrix formulation, it establishes desired variable constraints with good computational efficiency [11], as well as non-iterative design schemes based on linear matrix inequality (LMI) formulations, directly connected with stability and robustness [12]. Unfortunately, the majority of results related to those above-mentioned methods are not directly applicable to linear Metzler systems [13].

Because control algorithms used for Metzler linear positive systems are static in general, one of the motivating factors for the paper [14] was PID control law parameter design for single-input single-output (SISO) Metzler systems, but it turned out that its direct generalization for Metzler multiple-input multiple-output (MIMO) linear systems under given boundaries is not possible. To find out what underlying Metzler system restrictions must be placed on the system parameters, a particular structure of the matching conditions were analyzed and a generalized systematic scheme to determine the PID control for MIMO linear Metzlerian continuous-time systems is proposed in this paper. To show the Metzler and Hurwitz closed-loop system matrix structure, a more concrete PID control design task is formulated.

Concerning asymptotically stable solutions, the article is also a follow-up of another authors’ papers [15,16], which introduced a representation of Metzlerian system parameter constraints using LMIs, exploited a diagonal principle in stabilization and utilized a rhombic mapping of a strictly square Metzler matrix.

Such formulation is mathematically represented through the state-space models with a Metzler system matrix and supported by a minimal number of LMIs, defining structural constraints (compare [8]).

Reflecting the above-mentioned specific conditions, the approach presented in this paper for PID control design with application to linear Metzler MIMO systems is original and primary. Such inclusion of the D-part of PID, if the matching conditions are satisfied, refine system matrix parameter constraints in the design task. This is accomplished by assuming that a suitable equivalent system exists and the resulting PI design bilinear task can be tackled using LMIs and a linear matrix equality (LME) approach.

The remaining part of this paper is organized as follows. To present the reasoning path, brief comments on MIMO linear Metzlerian systems are given in Section 2. The proposed LMI technique, enforcing conditions on the PID control law design, with the structure of the matching conditions and the main theorem characterizing the system’s behavior, is stated in Section 3. More concretely, this section is focused on parametric features in PID control design for MIMO Metzlerian systems, basic constitutive relations concerning the D-part of the control law, feasibility problems that involve system parametric constraints and the ways of turning the approach into an LMI-based design formulation, conditioned by one LME. Confirming these results, Section 4 follows with an illustrative numerical example. Finally, Section 5 discusses the results and their interpretation to establish a straightforward perspective on the conclusions presented in Section 6.

Throughout the paper, the following notations are used: ${\mathit{x}}^{\mathrm{T}}$, ${\mathit{X}}^{\mathrm{T}}$ denotes the transpose of the vector x, and the matrix X, respectively, the indication ${\mathit{X}}^{h\mathrm{T}}$ means transpose of the h-th power of a square matrix X, the notation $\mathit{X}\otimes \mathit{Y}$ represents the Kronecker product (tensor product) of two real matrices $\mathit{X},\mathit{Y}$, $diag[·]$ outlines a diagonal form, $\rho \left(\mathit{X}\right)$ identifies a set of related matrix eigenvalues, labeling of matrix $\mathit{X}\succ 0$ means its positive definiteness, ${\mathit{I}}_n\in {\mathbb{R}}^{n\times n}$ is a unit matrix, $a\in {\mathbb{R}}_{+}$ is a non-negative real scalar, $\left({\mathbb{R}}_{+}^{n\times r}\right)$, ${\mathbb{R}}^{n\times r}$ refers to the set of $n\times r$ (non-negative) real matrices and ${\mathbb{M}}_{-+}^{n\times n}$, (${\mathbb{M}}_{-+\circ }^{n\times n}$) means the set of strictly (purely) Metzler square matrices, respectively.

2. Linear Metzler Systems Formalism and Control System Strategies

Making additional assumptions on the control design for a multiple-input, multiple-output (MIMO) Metzler positive system, the goal of this section is to present coincided conditions, which can be justified when defining this task as the problem of synthesis with the set of parametric constraints.

A linear, time-invariant continuous-time MIMO Metzler positive system allows the state-space description

$$\dot{\mathit{q}}\left(t\right)=\mathit{A}\mathit{q}\left(t\right)+\mathit{B}\mathit{u}\left(t\right),$$

(1)

$$\mathit{y}\left(t\right)=\mathit{C}\mathit{q}\left(t\right),$$

(2)

where $\mathit{q}\left(t\right)\in {\mathbb{R}}_{+}^n$, $\mathit{u}\left(t\right)\in {\mathbb{R}}^r$, $\mathit{y}\left(t\right)\in {\mathbb{R}}_{+}^m$ are the system state vector, control input and measurable output.

Since there exist different techniques to preserve important properties of the Metzler linear positive systems for parameter constraints compartmentalization, with reference to the notation, it is considered that $\mathit{B}\in {\mathbb{R}}_{+}^{n\times m}$, $\mathit{C}\in {\mathbb{R}}_{+}^{m\times n}$ are non-negative matrices and $\mathit{A}\in {\mathbb{M}}_{-+}^{n\times n}$ is Metzler.

Exploiting the Metzler matrix structure notation, the following features have to be highlighted.

Definition 1

([17]). A square matrix $\mathit{A}\in {\mathbb{M}}_{-+\circ }^{n\times n}$ is purely Metzler if its diagonal elements are negative and its off-diagonal elements are non-negative. A square matrix $\mathit{A}\in {\mathbb{M}}_{-+}^{n\times n}$ is strictly Metzler if its diagonal elements are negative and its off-diagonal elements are positive. A Metzler matrix is stable if it is Hurwitz. From a strictly Metzler matrix $\mathit{A}\in {\mathbb{M}}_{-+}^{n\times n}$ imply $n^2$ structural constraints

$$a_{ii}<0\forall i=1,\cdots n,a_{ij,i\ne j}>0\forall i,j=1,\cdots n.$$

(3)

Remark 1.

Since $\mathit{B}\in {\mathbb{R}}_{+}^{n\times r}$, $\mathit{C}\in {\mathbb{R}}_{+}^{m\times n}$ are non-negative, a negative feedback makes smaller (non-negative or positive) off-diagonal elements, and it could destroy the Metzler structure, setting one or more off-diagonal elements to a negative value. This fact also highlights that structural constraints must be included in the synthesis conditions to keep the resulting Metzler structure.

While for general linear systems it is possible to work with signum indefinite elements in the matrix inversion of a square matrix, for Metzler systems it may be difficult, or impossible, to provide general statements if this matrix operation has to be performed. Since a square matrix X and its inverse have non-negative, structure if X is positively definite diagonal, to guarantee structural constraints the LMI based design conditions for Metzler systems are formulated using positive definite diagonal matrix variables, and the term "diagonal stability" is used [15,18]. If $\mathit{A}\in {\mathbb{R}}_{-+}^{n\times n}$ is only purely Metzler, the synthesis conditions have to reflect further structural constraints, includable in the design by related structured diagonal matrix variables [19].

Proposition 1

([1]). A solution $\mathit{q}\left(t\right)$ of (1) for $t\ge 0$ is asymptotically stable and positive if $\mathit{A}\in {\mathbb{M}}_{-+}^{n\times n}$ is a stable Metzler matrix, $\mathit{B}\in {\mathbb{R}}_{+}^{n\times r}$ is a non-negative matrix and the state vector $\mathit{q}\left(t\right)\in {\mathbb{R}}_{+}^n$ for given $\mathit{u}\left(t\right)\in {\mathbb{R}}_{+}$ and $\mathit{q}\left(0\right)\in {\mathbb{R}}_{+}$. The linear system (1), (2) is asymptotically stable and positive if $\mathit{A}\in {\mathbb{M}}_{-+}^{n\times n}$ is a stable Metzler matrix, $\mathit{B}\in {\mathbb{R}}_{+}^{n\times r}$, $\mathit{C}\in {\mathbb{R}}_{+}^{m\times n}$ are non-negative matrices and the output vector $\mathit{y}\left(t\right)\in {\mathbb{R}}_{+}$ for all $\mathit{u}\left(t\right)\in {\mathbb{R}}_{+}$ and $\mathit{q}\left(0\right)\in {\mathbb{R}}_{+}$.

Definition 2

([20]). A matrix $\mathit{L}\in {\mathbb{R}}^{n\times n}$ is a permutation matrix if exactly one item in each column and row is equal to 1 and all other elements are equal to 0.

Taking into account Definition 2 and envisaging a diagonal $\mathit{Y}\in {\mathbb{R}}^{n\times n}$ such that

$$\mathit{Y}=\mathrm{diag}\left[\begin{array}{cccc}{y}_1& y_2& \cdots & y_n\end{array}\right],$$

(4)

then it yields

$${\mathit{L}}^{\mathrm{T}}\mathit{Y}\mathit{L}=\mathrm{diag}\left[\begin{array}{cccc}{y}_2& \cdots & y_n& y_1\end{array}\right],$$

(5)

if ${\mathit{L}}^{\mathrm{T}}\in {\mathbb{R}}^{n\times n}$ takes the circulant form

$${\mathit{L}}^{\mathrm{T}}=\left[\begin{array}{cc}0& {\mathit{I}}_{n-1}\\ 1& 0\end{array}\right].$$

(6)

Remark 2

([16]). The diagonal stabilization problem can be reformulated using a rhombic mapping of the square strictly Metzler matrix $\mathit{A}\in {\mathbb{M}}_{-+}^{n\times n}$ of the form

$$\mathit{A}=\left[\begin{array}{ccccc}{a}_{11}& a_{12}& a_{13}& \cdots & a_{1n}\\ a_{21}& a_{22}& a_{23}& \cdots & a_{2n}\\ a_{31}& a_{32}& a_{33}& \cdots & a_{3n}\\ ⋮& ⋮& ⋮& ⋮& ⋮\\ a_{n1}& a_{n2}& a_{n3}& \cdots & a_{nn}\end{array}\right],$$

(7)

where the rhombic mapping is constructed using circular shifts of rows of $\mathit{A}\in {\mathbb{M}}_{-+}^{n\times n}$ as

$${\mathit{A}}_{\mathsf{\Theta }}=\left[\begin{array}{cccccc}{a}_{11}& a_{12}& a_{13}& \cdots & a_{1n}& \\ & a_{22}& a_{23}& \cdots & a_{2n}& a_{21}& \\ & & a_{33}& \cdots & a_{3n}& a_{31}& a_{32}\\ & & & \ddots & ⋮& ⋮& ⋮& \ddots \\ & & & & a_{nn}& a_{n1}& a_{n2}& \cdots & a_{n,n-1}\end{array}\right].$$

(8)

It is evident that generally $n^2$ parametric constraints (3) can be defined by the negativeness of ${\mathit{A}}_{\mathsf{\Theta }}(i,i+h)\in {\mathbb{R}}_{+}^{n\times n}$ for $h=0$ and by the positiveness of $(n-1)$ diagonal matrices ${\mathit{A}}_{\mathsf{\Theta }}(i,i+h)\in {\mathbb{R}}_{+}^{n\times n}$ for $h=1,\cdots ,n-1$ with

$${\mathit{A}}_{\mathsf{\Theta }}(i,i+h)=\mathit{diag}\left[\begin{array}{cccccc}{a}_{1,1+h}& \cdots & a_{n-h,n}& a_{n-h+1,1}& \cdots & a_{n,h}\end{array}\right]$$

(9)

related to the diagonals of (8).

Definition 3

([17]). Let $\mathit{U}\in {\mathbb{R}}^{m\times m}$, $\mathit{V}\in {\mathbb{R}}^{n\times n}$ then the $\left(mn\right)$-dimensional matrix, called the Kronecker product of U and V, is constructed as

$$\mathit{U}\otimes \mathit{V}=\left[{\left\{u_{ij}\mathit{V}\right\}}_{i,j=1}^m\right],\mathit{U}=\left[{\left\{u_{ij}\right\}}_{i,j=1}^m\right].$$

(10)

It can be underlined at this point that the following Kronecker product properties [21] will now be priority

$$({\mathit{I}}_n\otimes \mathit{U})(\mathit{V}\otimes {\mathit{I}}_m)=(\mathit{V}\otimes {\mathit{I}}_m)({\mathit{I}}_n\otimes \mathit{U}),$$

(11)

$${(\mathit{U}\otimes \mathit{V})}^{-1}={\mathit{U}}^{-1}\otimes {\mathit{V}}^{-1},$$

(12)

$${(\mathit{U}\otimes \mathit{V})}^{\mathrm{T}}={\mathit{U}}^{\mathrm{T}}\otimes {\mathit{V}}^{\mathrm{T}}.$$

(13)

Consider the system (1), (2) and the properties of diagonals of the mapping (8), (9) with the specific relation to the Metzler system diagonal stabilization principle. To keep the notation simple, without loss of generality, this principle is briefly formulated in the following lemma.

Lemma 1

([22]). Let a square real $n\times n$ matrix Λ be partitioned as

$$\mathsf{\Lambda }=\mathit{A}-\mathit{B}\mathit{D}\mathit{C},$$

(14)

where $\mathit{A}\in {\mathbb{M}}_{-+}^{n\times n}$, $\mathit{B}\in {\mathbb{R}}_{+}^{n\times m}$, $\mathit{C}\in {\mathbb{R}}_{+}^{m\times n}$, $\mathit{D}\in {\mathbb{R}}_{+}^{m\times m}$, while A is strictly Metzler. Then Λ is strictly Metzler if, equivalently,

(i) 

$$\begin{array}{cc}\hfill a_{ii}-{\mathit{b}}_i^{\mathrm{T}}\mathit{D}{\mathit{c}}_i<0& \mathit{for}\mathit{all}i=1,\cdots ,n,\hfill \\ \hfill a_{ij}-{\mathit{b}}_i^{\mathrm{T}}\mathit{D}{\mathit{c}}_j>0& \mathit{for}\mathit{all}i,j=1,\cdots n,i\ne j,\hfill \end{array}$$

(15)

(ii) 

$$\begin{array}{cc}\hfill & {\mathit{A}}_{\mathsf{\Theta }}(i,i)-{\mathit{B}}_d{\mathit{D}}_d{\mathit{C}}_d\prec 0,\hfill \\ \hfill & {\mathit{A}}_{\mathsf{\Theta }}(i,i+h)-{\mathit{B}}_d{\mathit{D}}_d{\mathit{C}}_{dh}\succ 0,\hfill \end{array}$$

(16)

where

$${\mathit{B}}^{\mathrm{T}}=\left[\begin{array}{ccc}{\mathit{b}}_1& \cdots & {\mathit{b}}_n\end{array}\right],{\mathit{B}}_d=\mathit{diag}\left[\begin{array}{ccc}{\mathit{b}}_1^{\mathrm{T}}& \cdots & {\mathit{b}}_n^{\mathrm{T}}\end{array}\right],$$

(17)

$$\mathit{C}=\left[\begin{array}{ccc}{\mathit{c}}_1& \cdots & {\mathit{c}}_n\end{array}\right],{\mathit{C}}_d=\mathit{diag}\left[\begin{array}{ccc}{\mathit{c}}_1& \cdots & {\mathit{c}}_n\end{array}\right],$$

(18)

$${\mathit{D}}_d={\mathit{I}}_n\otimes \mathit{D},{\mathit{C}}_{dh}={\mathit{S}}^{h\mathrm{T}}{\mathit{C}}_d{\mathit{L}}^h,\mathit{S}=\mathit{L}\otimes {\mathit{I}}_m.$$

(19)

Moreover, the square matrix representation from its rhombic diagonals is given as

$$\mathsf{\Lambda }=\sum _{h=0}^{n-1}({\mathit{A}}_{\mathsf{\Theta }}(i,i+h)-{\mathit{B}}_d{\mathit{D}}_d{\mathit{C}}_{dh}){\mathit{L}}^{h\mathrm{T}}.$$

(20)

Since the positivity of the systems is defined by a non-negative system state, the non-negative system input and output matrix parameters and a Metzler system matrix structure, it is necessary to take these facts into account when synthesizing the PID MIMO controller.

Because this introduces an added limitation in the synthesis conditions, the aim is to develop a systematic framework for PID control design for a given class of Metzler linear MIMO systems that would be sufficiently general with respect to the system input matrix parameter and the related structure of the measured system but also effective in terms of closed-loop system stability and positivity. The latter is a multi-variable problem, subject to the given parametric constraints.

The problem is finally the following: Assuming the given class of Metzler linear MIMO systems, the matching conditions for the existence of PID control with transition to the input and output matrices and the design procedure of the Metzler system have to be defined, based on the general set of LMIs that is given if the matching conditions are satisfied. Since the use of other known PID control design formalism for this class of system remains rather heuristic in the system positiveness, the obtained results cannot be compared with unknown positive solutions.

3. Main Results

To provide a constructive solution to the parameter feasibility problem in the synthesis of PID controllers for a given system Metzler system class, it is necessary to establish the matching conditions related to the system input matrix and the system output matrix, the LMI representation of the Metzler system matrix parametric constraints and a direct consequence of the control law matrix parameters on the closed-loop system asymptotic stability and positiveness.

3.1. Parametric Features in PID Control Design

To respect the positiveness of the system variables for the considered class of square Metzler linear MIMO systems (1), (2), the MIMO continuous-time PID control algorithm is considered as

$$\mathit{u}\left(t\right)={\mathit{K}}_P{\mathit{e}}_r\left(t\right)+{\mathit{K}}_I{\int }_0^t{\mathit{C}}_p\mathit{q}\left(\tau \right)d\tau -{\mathit{K}}_D{\dot{\mathit{e}}}_r\left(t\right),$$

(21)

where $r=m$, $w_r\in {\mathbb{R}}_{+}^m$ is a constant positive reference output vector, ${\mathit{e}}_r\left(t\right)\in {\mathbb{R}}^m$ is the control signal error vector, where

$${\mathit{e}}_r\left(t\right)={\mathit{w}}_r-\mathit{y}\left(t\right)$$

(22)

and ${\mathit{K}}_P,{\mathit{K}}_I,{\mathit{K}}_D\in {\mathbb{R}}_{+}^{r\times m}$ are non-negative matrix parameters of the controller.

To reduce the desired specification on the closed-loop system, it is assumed that there exist associated $\mathit{w}\in {\mathbb{R}}_{+}^n$ and $\mathit{e}\left(t\right)\in {\mathbb{R}}^n$ such that

$${\mathit{w}}_r=\mathit{C}\mathit{w},\mathit{e}\left(t\right)=\mathit{w}-\mathit{q}\left(t\right),{\mathit{e}}_r\left(t\right)=\mathit{C}\mathit{e}\left(t\right)$$

(23)

and the essential feature of the considered synthesis problem depending (23) is given as

$$\begin{array}{cc}\hfill \mathit{u}\left(t\right)& ={\mathit{K}}_P\mathit{C}\mathit{e}\left(t\right)+{\mathit{K}}_I\mathit{p}\left(t\right)-{\mathit{K}}_D\mathit{C}\dot{\mathit{e}}\left(t\right)\hfill \\ \hfill & =-{\mathit{K}}_P\mathit{C}\mathit{q}\left(t\right)+{\mathit{K}}_I\mathit{p}\left(t\right)+{\mathit{K}}_D\mathit{C}\dot{\mathit{q}}\left(t\right)+{\mathit{K}}_P\mathit{C}\mathit{w}\hfill \\ \hfill & =-{\mathit{K}}_P\mathit{C}\mathit{q}\left(t\right)+{\mathit{K}}_I\mathit{p}\left(t\right)+{\mathit{K}}_D\mathit{C}\dot{\mathit{q}}\left(t\right)+{\mathit{K}}_P{\mathit{w}}_r,\hfill \end{array}$$

(24)

while the variable vector $\dot{\mathit{p}}\left(t\right)\in {\mathbb{R}}_{+}^m$ on input of the integrator is

$$\dot{\mathit{p}}\left(t\right)\left(t\right)={\mathit{C}}_p\mathit{q}\left(t\right).$$

(25)

Note, in general, it can set ${\mathit{C}}_p=\mathit{C}$ and, appending to the integrator, input all state variables involved in the measurable system output projection.

For actual computations, the equivalence of the assembled system structure can be expressed well as

$$\left[\begin{array}{cc}{\mathit{I}}_n-\mathit{B}{\mathit{K}}_D\mathit{C}& 0\\ 0& {\mathit{I}}_m\end{array}\right]\left[\begin{array}{c}\dot{\mathit{q}}\left(t\right)\\ \dot{\mathit{p}}\left(t\right)\end{array}\right]=\left[\begin{array}{cc}\mathit{A}-\mathit{B}{\mathit{K}}_P\mathit{C}& \mathit{B}{\mathit{K}}_I\\ {\mathit{C}}_p& -{\mathit{I}}_m\end{array}\right]\left[\begin{array}{c}\mathit{q}\left(t\right)\\ \mathit{p}\left(t\right)\end{array}\right]+\left[\begin{array}{cc}\mathit{B}{\mathit{K}}_P& 0\\ 0& {\mathit{I}}_m\end{array}\right]\left[\begin{array}{c}{\mathit{w}}_r\\ \mathit{p}\left(t\right)\end{array}\right],$$

(26)

$$\left[\begin{array}{c}\mathit{y}\left(t\right)\\ \mathit{p}\left(t\right)\end{array}\right]=\left[\begin{array}{cc}\mathit{C}& 0\\ 0& {\mathit{I}}_m\end{array}\right]\left[\begin{array}{c}\mathit{q}\left(t\right)\\ \mathit{p}\left(t\right)\end{array}\right].$$

(27)

In the above view, the composite variables and the composite matrix parameters are introduced as

$${\mathit{q}}^{\circ }\left(t\right)=\left[\begin{array}{c}\mathit{q}\left(t\right)\\ \mathit{p}\left(t\right)\end{array}\right],{\mathit{w}}^{\circ }\left(t\right)=\left[\begin{array}{c}{\mathit{w}}_r\\ \mathit{p}\left(t\right)\end{array}\right],{\mathit{y}}^{\circ }\left(t\right)=\left[\begin{array}{c}\mathit{y}\left(t\right)\\ \mathit{p}\left(t\right)\end{array}\right],$$

(28)

$${\mathit{A}}^{\diamond }=\left[\begin{array}{cc}\mathit{A}-\mathit{B}{\mathit{K}}_P\mathit{C}& \mathit{B}{\mathit{K}}_I\\ {\mathit{C}}_p& -{\mathit{I}}_m\end{array}\right],{\mathit{B}}^{\diamond }=\left[\begin{array}{cc}\mathit{B}{\mathit{K}}_P& 0\\ 0& {\mathit{I}}_m\end{array}\right],{\mathit{E}}^{\diamond }=\left[\begin{array}{cc}{\mathit{I}}_n-\mathit{B}{\mathit{K}}_D\mathit{C}& 0\\ 0& {\mathit{I}}_m\end{array}\right],{\mathit{C}}^{\circ }=\left[\begin{array}{cc}\mathit{C}& 0\\ 0& {\mathit{I}}_m\end{array}\right],$$

(29)

which specify the corresponding closed loop system description

$${\mathit{E}}^{\diamond }{\dot{\mathit{q}}}^{\circ }\left(t\right)={\mathit{A}}^{\diamond }{\mathit{q}}^{\circ }\left(t\right)+{\mathit{B}}^{\diamond }{\mathit{w}}^{\circ }\left(t\right),$$

(30)

$${\mathit{y}}^{\circ }\left(t\right)={\mathit{C}}^{\circ }{\mathit{q}}^{\circ }\left(t\right).$$

(31)

Thus, if ${\mathit{E}}^{\diamond }$ is regular, then

$${\dot{\mathit{q}}}^{\circ }\left(t\right)={\mathit{E}}^{\diamond -1}{\mathit{A}}^{\diamond }{\mathit{q}}^{\circ }\left(t\right)+{\mathit{E}}^{\diamond -1}{\mathit{B}}^{\diamond }{\mathit{w}}^{\circ }\left(t\right).$$

(32)

With this non-descriptor notation, the matrix inequality procedures can hold true but the direct convexifying lead to bi-linear matrix inequalities.

3.2. Basic Constitutive Control Constraints

To search for a stabilizing PID MIMO controller in which the design is LMI-able, the role of the system matrix parameter in the design is analyzed.

If $\mathit{A}\in {\mathbb{M}}_{-+}^{n\times n}$, $\mathit{B}\in {\mathbb{R}}_{+}^{n\times r}$, $\mathit{C}\in {\mathbb{R}}_{+}^{m\times n}$ are the parameters of a Metzler positive MIMO system, the matrices ${\mathit{B}}^{\diamond }$, ${\mathit{C}}^{\circ }$ are non-negative. Then, with non-negative parameters of the PID controller ${\mathit{K}}_P,{\mathit{K}}_I,{\mathit{K}}_D\in {\mathbb{R}}_{+}^{r\times m}$ the matrix ${\mathit{A}}^{\diamond }$ is Metzler if $(\mathit{A}-\mathit{B}{\mathit{K}}_P{\mathit{C}}^{\mathrm{T}})$ is (strictly) Metzler, $({\mathit{I}}_n-\mathit{B}{\mathit{K}}_D\mathit{C})$ is regular and ${({\mathit{I}}_n-\mathit{B}{\mathit{K}}_D\mathit{C})}^{-1}$ is positive.

The following simplification is used to show how to realize positive effects of the PID controller derivative part on the Metzler structural constraints satisfying the above assumption.

Consider a sub-class of square linear Metzler MIMO systems, where $r=m$, $\mathit{A}\in {\mathbb{M}}_{-+}^{n\times n}$ is strictly Metzler, and $\mathit{B}\in {\mathbb{R}}_{+}^{n\times m}$, $\mathit{C}\in {\mathbb{R}}_{+}^{m\times n}$ take the following non-negative structures

$$\mathit{C}=\left[\begin{array}{cc}{\mathit{I}}_m& 0\end{array}\right],{\mathit{B}}^{\mathrm{T}}=\left[\begin{array}{cc}{\mathit{B}}_1^{\mathrm{T}}& {\mathit{B}}_2^{\mathrm{T}}\end{array}\right],$$

(33)

while ${\mathit{B}}_1\in {\mathbb{R}}_{+}^{m\times m}$ is regular and non-negative. If ${\mathit{K}}_D\in {\mathbb{R}}_{+}^{m\times m}$ is supposed as regular and non-negative, using the Sherman–Morrison–Woodbury formula [23] yields

$$\begin{array}{cc}\hfill {({\mathit{I}}_n-\mathit{B}{\mathit{K}}_D\mathit{C})}^{-1}& ={\mathit{I}}_n-\mathit{B}{({\mathit{K}}_D^{-1}+\mathit{C}\mathit{B})}^{-1}\mathit{C}\hfill \\ \hfill & ={\mathit{I}}_n-\mathit{B}{(-{\mathit{K}}_D^{-1}+{\mathit{B}}_1)}^{-1}\mathit{C}.\hfill \end{array}$$

(34)

Setting with a positive scalar $\epsilon \in {\mathbb{R}}_{+}$ that

$${\mathit{K}}_D^{-1}=(1+{\epsilon }^{-1}){\mathit{B}}_1,$$

(35)

then

$${({\mathit{I}}_n-\mathit{B}{\mathit{K}}_D\mathit{C})}^{-1}={\mathit{I}}_n-\mathit{B}{(-{\epsilon }^{-1}{\mathit{B}}_1)}^{-1}\mathit{C}={\mathit{I}}_n+\epsilon \mathit{B}{\mathit{B}}_1^{-1}\mathit{C}$$

(36)

and the block structure of B implies

$${({\mathit{I}}_n-\mathit{B}{\mathit{K}}_D\mathit{C})}^{-1}={\mathit{I}}_n+\epsilon \left[\begin{array}{c}{\mathit{B}}_1\\ {\mathit{B}}_2\end{array}\right]{\mathit{B}}_1^{-1}\mathit{C}={\mathit{I}}_n+\epsilon \left[\begin{array}{c}\mathit{C}\\ {\mathit{B}}_2{\mathit{B}}_1^{-1}\mathit{C}\end{array}\right].$$

(37)

Since the inverse of a square positive matrix is, in general, signum indefinite, it is evident that ${\mathit{B}}_1$ must have a non-negative structure, at least such that one of its non-diagonal elements takes the value of zero, or ${\mathit{B}}_1$ has to be a diagonal positive definite matrix. Under these conditions, ${\mathit{B}}_1^{-1}$ is non-negative and it yields

$${({\mathit{I}}_n-\mathit{B}{\mathit{K}}_D\mathit{C})}^{-1}={\mathit{I}}_n+\epsilon \mathit{B}{\mathit{B}}_1^{-1}\mathit{C}\succ 0,$$

(38)

$${\mathit{E}}^{\diamond -1}={\left[\begin{array}{cc}{\mathit{I}}_n-\mathit{B}{\mathit{K}}_D\mathit{C}& 0\\ 0& {\mathit{I}}_m\end{array}\right]}^{-1}=\left[\begin{array}{cc}{\mathit{I}}_n+\epsilon \mathit{B}{\mathit{B}}_1^{-1}\mathit{C}& 0\\ 0& {\mathit{I}}_m\end{array}\right]\succ 0,$$

(39)

respectively.

Remark 3.

If B is positive, with a non-negative C of full rank, only the positive ${\mathit{B}}_1$ can be constructed and the positive ${\mathit{E}}^{\diamond -1}$ cannot be obtained. The case when the matrix B is non-negative, ${\mathit{B}}_1$ is positive and ${\mathit{B}}_2=0$ gives the positive ${\mathit{E}}^{\diamond -1}$ but a signum indefinite control parameter ${\mathit{K}}_D$.

The positive effect of the PID controller derivative part on the Metzler matrix structure can be exploited if C is of full rank and non-negative, the matrix B is non-negative and ${\mathit{B}}_1$ is diagonal.

An ad hoc solution can exist if ${({\mathit{I}}_n-\mathit{B}{\mathit{K}}_D\mathit{C})}^{-1}$ is not positive and the structure of ${\mathit{B}}_1$ is, e.g., a degenerative lower triangular matrix, but it hardly depends on the parameters of Metzler matrix A.

The main idea of this remark can be generalized with respect to the diagonal stabilization principle by the following matching conditions.

Definition 4.

The derivative part of the PID controller in control of strictly Metzler linear MIMO systems exists if non-negative matrices $\mathit{B}\in {\mathbb{R}}_{+}^{n\times m}$, $\mathit{C}\in {\mathbb{R}}_{+}^{m\times n}$ of the system satisfy the matching condition

$$\mathit{C}\mathit{B}={\mathit{B}}_1=diag\left[\begin{array}{cccc}{\mathit{b}}_{11}& {\mathit{b}}_{22}& \cdots & {\mathit{b}}_{mm}\end{array}\right]\succ 0,$$

(40)

where $\mathit{C}\mathit{B}\in {\mathbb{R}}_{+}^{m\times m}$, $\mathit{C}=\left[{\mathit{I}}_{m}0\right]$.

The matrix ${\mathit{B}}_1$, if exists, can be constructed by a linear coordinate transform of the system state variables while the associated measured variables must be chosen in such a way that C takes the above structure.

Remark 4.

Prescribing (35) as a solution to the positiveness problem, the direct evaluation then implies the diagonal control law parameter

$${\mathit{K}}_D=\frac{\epsilon }{1+\epsilon }{\mathit{B}}_1^{-1},$$

(41)

with dependence on a positive tuning parameter $\epsilon \in {\mathbb{R}}_{+}$ and positive definite diagonal ${\mathit{B}}_1$.

The problem can be put in the matrix form when the availing structure of the static output control is readily solved by the parametrization

$$\begin{array}{cc}\hfill {\mathit{A}}^{\diamond }& =\left[\begin{array}{cc}\mathit{A}& 0\\ {\mathit{C}}_p& -{\mathit{I}}_m\end{array}\right]-\left[\begin{array}{cc}\mathit{B}& 0\\ 0& {\mathit{I}}_m\end{array}\right]\left[\begin{array}{cc}{\mathit{K}}_P& -{\mathit{K}}_I\\ 0& 0\end{array}\right]\left[\begin{array}{cc}\mathit{C}& 0\\ 0& {\mathit{I}}_m\end{array}\right]\hfill \\ \hfill & ={\mathit{A}}^{\circ }-{\mathit{B}}^{\circ }{\mathit{K}}^{\circ }{\mathit{C}}^{\circ }\hfill \end{array}$$

(42)

where

$${\mathit{A}}^{\circ }=\left[\begin{array}{cc}\mathit{A}& 0\\ {\mathit{C}}_p& -{\mathit{I}}_m\end{array}\right],{\mathit{B}}^{\circ }=\left[\begin{array}{cc}\mathit{B}& 0\\ 0& {\mathit{I}}_m\end{array}\right],{\mathit{C}}^{\circ }=\left[\begin{array}{cc}\mathit{C}& 0\\ 0& {\mathit{I}}_m\end{array}\right],{\mathit{K}}^{\circ }=\left[\begin{array}{cc}{\mathit{K}}_P& -{\mathit{K}}_I\\ 0& 0\end{array}\right].$$

(43)

Therefore,

$${\mathit{E}}^{\diamond -1}{\mathit{B}}^{\diamond }=\left[\begin{array}{cc}{\mathit{I}}_n+\epsilon \mathit{B}{\mathit{B}}_1^{-1}\mathit{C}& 0\\ 0& {\mathit{I}}_m\end{array}\right]\left[\begin{array}{cc}\mathit{B}{\mathit{K}}_P& 0\\ 0& {\mathit{I}}_m\end{array}\right]=\left[\begin{array}{cc}\mathit{B}{\mathit{K}}_P+\epsilon \mathit{B}{\mathit{B}}_1^{-1}\mathit{C}\mathit{B}{\mathit{K}}_P& 0\\ 0& {\mathit{I}}_m\end{array}\right]={\mathit{B}}_w^{•},$$

(44)

$${\mathit{E}}^{\diamond -1}{\mathit{B}}^{\circ }=\left[\begin{array}{cc}{\mathit{I}}_n+\epsilon \mathit{B}{\mathit{B}}_1^{-1}\mathit{C}& 0\\ 0& {\mathit{I}}_m\end{array}\right]\left[\begin{array}{cc}\mathit{B}& 0\\ 0& {\mathit{I}}_m\end{array}\right]=\left[\begin{array}{cc}\mathit{B}+\epsilon \mathit{B}{\mathit{B}}_1^{-1}\mathit{C}\mathit{B}& 0\\ 0& {\mathit{I}}_m\end{array}\right]={\mathit{B}}^{•},$$

(45)

$${\mathit{E}}^{\diamond -1}{\mathit{A}}^{\circ }=\left[\begin{array}{cc}{\mathit{I}}_n+\epsilon \mathit{B}{\mathit{B}}_1^{-1}\mathit{C}& 0\\ 0& {\mathit{I}}_m\end{array}\right]\left[\begin{array}{cc}\mathit{A}& 0\\ {\mathit{C}}_p& -{\mathit{I}}_m\end{array}\right]=\left[\begin{array}{cc}\mathit{A}+\epsilon \mathit{B}{\mathit{B}}_1^{-1}\mathit{C}\mathit{A}& 0\\ {\mathit{C}}_p& -{\mathit{I}}_m\end{array}\right]={\mathit{A}}^{•}.$$

(46)

It is evident, since B is non-negative that ${\mathit{B}}^{•}$ as well as ${\mathit{B}}_v^{•}$ are non-negative, and since A is strictly Metzler, then ${\mathit{A}}^{•}$ is purely Metzler.

Using the above, (32) can be transformed into the form

$${\dot{\mathit{q}}}^{\circ }\left(t\right)=({\mathit{A}}^{•}-{\mathit{B}}^{•}{\mathit{K}}^{\circ }{\mathit{C}}^{\circ }){\mathit{q}}^{\circ }\left(t\right)+{\mathit{B}}_w^{•}{\mathit{w}}^{\circ }\left(t\right)={\mathit{A}}_c^{•}{\mathit{q}}^{\circ }\left(t\right)+{\mathit{B}}_w^{•}{\mathit{w}}^{\circ }\left(t\right),$$

(47)

also emphasizing that such system description leads directly to a bilinear structure of matrix inequalities.

Remark 5.

Since the resulting ${\mathit{B}}^{•}$, ${\mathit{B}}_w^{•}$ are non-negative, this in turn means that for a positive diagonal ${\mathit{K}}_D$ defined as in (41) it is sufficient to include in the synthesis the parametric constraints resulting from the desired Metzler structure of ${\mathit{A}}_c^{•}={\mathit{A}}^{•}-{\mathit{B}}^{•}{\mathit{K}}_{PI}^{\circ }{\mathit{C}}^{\circ }$, where ${\mathit{A}}^{•}$ is purely Metzler.

In order to provide constraint limitations, constraint structures need to reflect linear matrix inequality forms, but the structure of ${\mathit{A}}_c^{•}$ implies an essentially bilinear matrix inequality formulation. This can be eliminated by applying one linear matrix equality into the design (see, e.g., [22]). The principle is explicated at the point of application.

3.3. PID Control Law Parameter Design

Exploiting the diagonal stabilization principle in accession to the control design for strictly linear Metzlerian structures [15], the following matrix parameter ${\mathit{A}}^{•}\in {\mathbb{M}}_{-+}^{(n+m)\times (n+m)}$, ${\mathit{B}}^{•}\in {\mathbb{R}}_{+}^{(n+m)\times 2m}$, ${\mathit{C}}^{\circ }\in {\mathbb{R}}_{+}^{2m\times (n+m)}$ needs to be adequately represented diagonally. These are represented in accordance with Lemma 1 as follows:

$${\mathit{A}}^{•}=\left[\begin{array}{ccccccc}\hfill -a_{11}^{•}& a_{12}^{•}& \cdots \hfill & a_{1n}^{•}& 0& \cdots & 0\\ \hfill a_{21}^{•}& -a_{22}^{•}& \cdots \hfill & a_{2n}^{•}& 0& \cdots & 0\\ \hfill ⋮& & \ddots \hfill & & \\ \hfill a_{n1}^{•}& a_{n2}^{•}& \cdots \hfill & -a_{nn}^{•}& 0& \cdots & 0\\ \hfill c_{p11}& c_{p12}& \cdots \hfill & c_{p1n}& -1& \cdots & 0\\ \hfill ⋮& & ⋮\hfill & & & \ddots & \\ \hfill c_{pr1}& c_{pr2}& \cdots \hfill & c_{prn}& 0& \cdots & -1\end{array}\right]=\left[\begin{array}{c}{\left\{a_{ij}^{•}\right\}}_{i,j=1}^{n+m}\end{array}\right],$$

(48)

$${\mathit{A}}_{\mathsf{\Theta }}^{•}(i,i)=\mathrm{diag}\left[\begin{array}{cccccc}-a_{1,1}^{•}& \cdots & -a_{nn}^{•}& -1& \cdots & -1\end{array}\right],$$

(49)

$${\mathit{A}}_{\mathsf{\Theta }}^{•}(i,i+h)=\mathrm{diag}\left[\begin{array}{cccccc}{a}_{1,1+h}^{•}& \cdots & a_{n+m-h,n+m}^{•}& a_{n+m-h+1,1}^{•}& \cdots & a_{n+m,h}^{•}\end{array}\right],h=1,\cdots ,n+m-1,$$

(50)

$${\mathit{B}}^{•}=\left[\begin{array}{cccccc}{\mathit{b}}_{11}^{•}& \cdots & {\mathit{b}}_{1m}^{•}& 0& \cdots & 0\\ ⋮& & ⋮& & \\ {\mathit{b}}_{n1}^{•}& \cdots & {\mathit{b}}_{nm}^{•}& 0& \cdots & 0\\ 0& \cdots & 0& 1& \cdots & 0\\ ⋮& & ⋮& & \ddots & \\ 0& \cdots & 0& 0& \cdots & 1\end{array}\right]=\left[\begin{array}{c}{\mathit{b}}_1^{•\mathrm{T}}\\ ⋮\\ {\mathit{b}}_n^{•\mathrm{T}}\\ {\mathit{b}}_{n+1}^{•\mathrm{T}}\\ ⋮\\ {\mathit{b}}_{n+m}^{•\mathrm{T}}\end{array}\right],{\mathit{B}}_d^{•}=\mathrm{diag}\left[\begin{array}{cccccc}{\mathit{b}}_1^{•\mathrm{T}}& \cdots & {\mathit{b}}_n^{•\mathrm{T}}& {\mathit{b}}_{n+1}^{•\mathrm{T}}& \cdots & {\mathit{b}}_{n+m}^{•\mathrm{T}}\end{array}\right],$$

(51)

$${\mathit{C}}^{\circ }=\left[\begin{array}{cc}\mathit{C}& 0\\ 0& {\mathit{I}}_m\end{array}\right]=\left[\begin{array}{cccccc}{\mathit{c}}_1^{\circ }& \cdots & {\mathit{c}}_n^{\circ }& {\mathit{c}}_{n+1}^{\circ }& \cdots & {\mathit{c}}_{n+m}^{\circ }\end{array}\right],{\mathit{C}}_d^{•}=\mathrm{diag}\left[\begin{array}{cccccc}{\mathit{c}}_1^{\circ }& \cdots & {\mathit{c}}_n^{\circ }& {\mathit{c}}_{n+1}^{\circ }& \cdots & {\mathit{c}}_{n+m}^{\circ }\end{array}\right].$$

(52)

To re-make the design construction in linear matrix inequality forms from the given bilinear structure, the related matrices are constructed in this way

$$\mathit{L}=\left[\begin{array}{cc}{0}^{\mathrm{T}}& 1\\ {\mathit{I}}_{n+m-1}& 0\end{array}\right],\mathit{S}=\mathit{L}\otimes {\mathit{I}}_{2m},{\mathit{J}}^{\mathrm{T}}=\left[\begin{array}{ccc}{\mathit{I}}_{2m}& \cdots & {\mathit{I}}_{2m}\end{array}\right],{\mathit{C}}_{dh}^{•}={\mathit{S}}^{h\mathrm{T}}{\mathit{C}}_d^{•}{\mathit{L}}^h$$

(53)

and the associated block-diagonal gain representation matrix

$${\mathit{K}}^{\circ }=\left[\begin{array}{cc}{\mathit{K}}_P& -{\mathit{K}}_I\\ 0& 0\end{array}\right],{\mathit{K}}_d^{•}=\mathrm{diag}\left[\begin{array}{ccc}{\mathit{K}}^{\circ }& \cdots & {\mathit{K}}^{\circ }\end{array}\right]={\mathit{I}}_{n+m}\otimes {\mathit{K}}^{\circ }.$$

(54)

is constructed via the same recipe.

Thus, the main result can now be presented and proven.

Theorem 1.

Let $\mathit{C}\in {\mathbb{R}}_{+}^{m\times n}$, $\mathit{B}\in {\mathbb{R}}_{+}^{n\times r}$ take the forms (33), where ${\mathit{B}}_1\in {\mathbb{R}}_{+}^{r\times r}$ is a regular positive definite diagonal matrix and for a given positive scalar $\epsilon \in {\mathbb{R}}_{+}$ the matrix ${\mathit{A}}^{•}\in {\mathbb{M}}_{-+}^{(n+m)\times (n+m)}$ is Metzler and ${\mathit{B}}^{•},{\mathit{B}}_w^{•}\in {\mathbb{R}}_{+}^{(n+m)\times 2m}$ are non-negative and ${\mathit{C}}^{\circ }\in {\mathbb{R}}_{+}^{2m\times (n+m)}$ is non-negative, then the closed-loop built on (1), (2) under PID control (21) is stable if there exist positive definite diagonal matrices $\mathit{P}\in {\mathbb{R}}_{+}^{(n+m)\times (n+m)}$, $\mathit{H}\in {\mathbb{R}}_{+}^{2m\times 2m}$ and a non-negative matrix $\mathit{R}\in {\mathbb{R}}_{+}^{2m\times 2m}$ such that for $h=1,2,\cdots n+m-1$,

$$\mathit{P}={\mathit{P}}^{\mathrm{T}}\succ 0,\mathit{H}={\mathit{H}}^{\mathrm{T}}\succ 0,$$

(55)

$${\mathit{A}}_{\mathsf{\Theta }}(i,i)\mathit{P}-{\mathit{B}}_d^{•}{\mathit{R}}_d^{•}{\mathit{C}}_d^{•}\prec 0,$$

(56)

$${\mathit{L}}^h{\mathit{A}}_{\mathsf{\Theta }}(i,i+h){\mathit{L}}^{h\mathrm{T}}\mathit{P}-{\mathit{L}}^h{\mathit{B}}_d^{•}{\mathit{S}}^{h\mathrm{T}}{\mathit{R}}_d^{•}{\mathit{C}}_d^{•}\succ 0,$$

(57)

$${\mathit{A}}^{•}\mathit{P}+\mathit{P}{\mathit{A}}^{•\mathrm{T}}-{\mathit{B}}_d^{•}{\mathit{R}}_d^{•}\mathit{J}{\mathit{J}}^{\mathrm{T}}{\mathit{C}}_d^{•}-{\mathit{C}}_d^{•\mathrm{T}}\mathit{J}{\mathit{J}}^{\mathrm{T}}{\mathit{R}}_d^{•\mathrm{T}}{\mathit{B}}_d^{•\mathrm{T}}\prec 0,$$

(58)

$${\mathit{C}}_d^{•}\mathit{P}={\mathit{H}}_d{\mathit{C}}_d^{•},$$

(59)

where

$$\mathit{R},{\mathit{R}}_d={\mathit{I}}_{n+m}\otimes \mathit{R},{\mathit{H}}_d={\mathit{I}}_{n+m}\otimes \mathit{H}.$$

(60)

are structured matrix variables.

Within a feasible solution with a suitable positive ε fixing a diagonal positive definite ${\mathit{K}}_D$, the gain ${\mathit{K}}^{\circ }\in {\mathbb{R}}_{+}^{2m\times 2m}$ representing not fixed design parameters is

$${\mathit{K}}^{\circ }=\left[\begin{array}{cc}{\mathit{K}}_P& -{\mathit{K}}_I\\ 0& 0\end{array}\right]=\mathit{R}{\mathit{H}}^{-1}.$$

(61)

Proof of Theorem 1.

For a stable realization of ${\mathit{A}}_c^{•}$, it yields according to Lyapunov inequality and relation (19)

$$\begin{array}{cc}\hfill & {\mathit{A}}_c^{•}\mathit{P}+\mathit{P}{\mathit{A}}_c^{•\mathrm{T}}=\hfill \\ \hfill & =\sum _{h=0}^n({\mathit{A}}_{\mathsf{\Theta }}^{•}(i,i+h){\mathit{L}}^{h\mathrm{T}}-{\mathit{B}}_d^{•}{\mathit{K}}_d^{•}{\mathit{C}}_{dh}^{•}{\mathit{L}}^{h\mathrm{T}})\mathit{P}+\sum _{h=0}^n\mathit{P}{({\mathit{A}}_{\mathsf{\Theta }}^{•}(i,i+h){\mathit{L}}^{h\mathrm{T}}-{\mathit{B}}_d^{•}{\mathit{K}}_d^{•}{\mathit{C}}_{dh}^{•}{\mathit{L}}^{h\mathrm{T}})}^{\mathrm{T}}\prec 0.\hfill \end{array}$$

(62)

Then, using the above (11) it follows that

$$\begin{array}{cc}\hfill {\mathit{B}}_d^{•}{\mathit{K}}_d^{•}{\mathit{C}}_{dh}^{•}{\mathit{L}}^{h\mathrm{T}}& ={\mathit{B}}_d^{•}{\mathit{K}}_d^{•}{\mathit{S}}^{h\mathrm{T}}{\mathit{C}}_d^{•}{\mathit{L}}^h{\mathit{L}}^{h\mathrm{T}}\hfill \\ \hfill & ={\mathit{B}}_d^{•}({\mathit{I}}_{n+m}\otimes {\mathit{K}}^{\circ })({\mathit{L}}^{h\mathrm{T}}\otimes {\mathit{I}}_{2m}){\mathit{C}}_d^{•}\hfill \\ \hfill & ={\mathit{B}}_d^{•}({\mathit{L}}^{h\mathrm{T}}\otimes {\mathit{I}}_{2m})({\mathit{I}}_{n+m}\otimes {\mathit{K}}^{\circ }){\mathit{C}}_d^{•}\hfill \\ \hfill & ={\mathit{B}}_d^{•}{\mathit{S}}^{h\mathrm{T}}{\mathit{K}}_d^{•}{\mathit{C}}_d^{•}.\hfill \end{array}$$

(63)

Since there is an admissible change of parameters, the product ${\mathit{K}}_d^{•}{\mathit{C}}_d^{•}$ can be written as follows

$$\begin{array}{cc}\hfill {\mathit{K}}_d^{•}{\mathit{C}}_d^{•}& =\left[\begin{array}{ccc}{\mathit{K}}^{\circ }\mathit{H}& & \\ & \ddots & \\ & & {\mathit{K}}^{\circ }\mathit{H}\end{array}\right]\left[\begin{array}{ccc}{\mathit{H}}^{-1}& & \\ & \ddots & \\ & & {\mathit{H}}^{-1}\end{array}\right]{\mathit{C}}_d^{•}\hfill \\ \hfill & ={\mathit{R}}_d{\mathit{H}}_d^{-1}{\mathit{C}}_d^{•},\hfill \end{array}$$

(64)

where

$$\mathit{R}={\mathit{K}}^{\circ }\mathit{H},{\mathit{R}}_d={\mathit{I}}_{n+m}\otimes {\mathit{K}}^{\circ }.$$

(65)

and making this substitution, (16) allows

$${\mathit{A}}_{\mathsf{\Theta }}^{•}(i,i)\mathit{P}-{\mathit{B}}_d^{•}{\mathit{K}}_d^{•}{\mathit{H}}_d{\mathit{H}}_d^{-1}{\mathit{C}}_d^{•}\mathit{P}\prec 0,$$

(66)

$${\mathit{A}}_{\mathsf{\Theta }}^{•}(i,i+h){\mathit{L}}^{hT}\mathit{P}-{\mathit{B}}_d^{•}{\mathit{S}}^{h\mathrm{T}}{\mathit{K}}_d^{•}{\mathit{H}}_d{\mathit{H}}_d^{-1}{\mathit{C}}_d^{•}\mathit{P}\succ 0.$$

(67)

It can be easily derived when prescribing

$${\mathit{H}}_d^{-1}{\mathit{C}}_d^{•}={\mathit{C}}_d^{•}{\mathit{P}}^{-1}$$

(68)

that (66), (67) imply (56), (57) and (68) gives (59), while the left multiplication of (67) by ${\mathit{L}}^h$ retains diagonal structures of LMIs.

To avoid additional structured variable’s phenomena in the design conditions [19], it can be taken as

$$\begin{array}{cc}\hfill {\mathit{A}}_c^{•}\mathit{P}+\mathit{P}{\mathit{A}}_c^{•\mathrm{T}}& =({\mathit{A}}^{•}-{\mathit{B}}^{•}{\mathit{K}}^{\circ }{\mathit{C}}^{\circ })\mathit{P}+\mathit{P}{({\mathit{A}}^{•}-{\mathit{B}}^{•}{\mathit{K}}^{\circ }{\mathit{C}}^{\circ })}^{\mathrm{T}}\hfill \\ \hfill & ={\mathit{A}}^{•}\mathit{P}+\mathit{P}{\mathit{A}}^{•\mathrm{T}}-{\mathit{B}}^{•}{\mathit{K}}^{\circ }{\mathit{C}}^{\circ }\mathit{P}-\mathit{P}{\mathit{C}}^{\circ \mathrm{T}}{\mathit{K}}^{\circ \mathrm{T}}{\mathit{B}}^{•\mathrm{T}}\hfill \\ \hfill & ={\mathit{A}}^{•}\mathit{P}+\mathit{P}{\mathit{A}}^{•\mathrm{T}}-{\mathit{B}}_d^{•}{\mathit{R}}_d\mathit{J}{\mathit{J}}^{\mathrm{T}}{\mathit{C}}_d^{•}-{\mathit{C}}_d^{•\mathrm{T}}\mathit{J}{\mathit{J}}^{\mathrm{T}}{\mathit{R}}_d^{\mathrm{T}}{\mathit{B}}_d^{•\mathrm{T}}\hfill \\ \hfill & \prec 0,\hfill \end{array}$$

(69)

Thus, (58) follows from the inequality (69) and verifies the system stability. This ends the proof. □

4. Illustrative Numerical Example

An Unstable Metzler linear MIMO system (1), (2) is considered to illustrate the design step with relation to the system matrix parameters

$$\mathit{A}=\left[\begin{array}{cccc}\hfill -3.380& \hfill 2.208& \hfill 4.715& \hfill 2.676\\ \hfill 1.881& \hfill -4.290& \hfill 2.050& \hfill 0.675\\ \hfill 2.067& \hfill 4.273& \hfill -6.654& \hfill 2.893\\ \hfill 1.148& \hfill 2.273& \hfill 1.343& \hfill -2.104\end{array}\right],\mathit{b}=\left[\begin{array}{cc}0.0410& 0\\ 0& 0.0203\\ 0.0114& 0.0315\\ 0.0114& 0.0170\end{array}\right],$$

$$\mathit{C}={\mathit{C}}_p=\left[\begin{array}{cccc}1& 0& 0& 0\\ 0& 1& 0& 0\end{array}\right],\mathit{C}\mathit{B}=\left[\begin{array}{cc}0.0410& 0\\ 0& 0.0203\end{array}\right].$$

Fixing the tuning parameter $\epsilon =0.01$ while solving the design problem results in

$${\mathit{K}}_D=\left[\begin{array}{cc}0.2415& 0\\ 0& 0.4877\end{array}\right],$$

$${\mathit{E}}^{\circ -1}=\left[\begin{array}{cccccc}1.0100& 0& 0& 0& 0& 0\\ 0& 1.0100& 0& 0& 0& 0\\ 0.0028& 0.0155& 1& 0& 0& 0\\ 0.0028& 0.0084& 0& 1& 0& 0\\ 0& 0& 0& 0& 1& 0\\ 0& 0& 0& 0& 0& 1\end{array}\right],{\mathit{A}}^{•}=\left[\begin{array}{cccccc}\hfill -3.4138& 2.2301& 4.7622\hfill & 2.7028& 0& 0\\ \hfill 1.8998& -4.3329& 2.0705\hfill & 0.6818& 0& 0\\ \hfill 2.0868& 4.2126& -6.6092\hfill & 2.9109& 0& 0\\ \hfill 1.1544& 2.2432& 1.3732\hfill & -2.0909& 0& 0\\ \hfill 1& 0& 0\hfill & 0& -1& 0\\ \hfill 0& 1& 0\hfill & 0& 0& -1\end{array}\right],$$

$${\mathit{B}}^{•}=\left[\begin{array}{cccc}0.0414& 0& 0& 0\\ 0& 0.0205& 0& 0\\ 0.0115& 0.0318& 0& 0\\ 0.0115& 0.0172& 0& 0\\ 0& 0& 1& 0\\ 0& 0& 0& 1\end{array}\right]=\left[\begin{array}{c}{\mathit{b}}_1^{•\mathrm{T}}\\ {\mathit{b}}_2^{•\mathrm{T}}\\ {\mathit{b}}_3^{•\mathrm{T}}\\ {\mathit{b}}_4^{•\mathrm{T}}\\ {\mathit{b}}_5^{•\mathrm{T}}\\ {\mathit{b}}_6^{•\mathrm{T}}\end{array}\right],{\mathit{B}}_d^{•}=\mathrm{diag}\left[\begin{array}{cccccc}{\mathit{b}}_1^{•\mathrm{T}}& {\mathit{b}}_2^{•\mathrm{T}}& {\mathit{b}}_3^{•\mathrm{T}}& {\mathit{b}}_4^{•\mathrm{T}}& {\mathit{b}}_5^{•\mathrm{T}}& {\mathit{b}}_6^{•\mathrm{T}}\end{array}\right],$$

$${\mathit{C}}^{\circ }=\left[\begin{array}{cccccc}1& 0& 0& 0& 0& 0\\ 0& 1& 0& 0& 0& 0\\ 0& 0& 0& 0& 1& 0\\ 0& 0& 0& 0& 0& 1\end{array}\right]=\left[\begin{array}{cccccc}{\mathit{c}}_1^{\circ }& {\mathit{c}}_2^{\circ }& {\mathit{c}}_3^{\circ }& {\mathit{c}}_4^{\circ }& {\mathit{c}}_5^{\circ }& {\mathit{c}}_6^{\circ }\end{array}\right],{\mathit{C}}_d^{•}=\mathrm{diag}\left[\begin{array}{cccccc}{\mathit{c}}_1^{\circ }& {\mathit{c}}_2^{\circ }& {\mathit{c}}_3^{\circ }& {\mathit{c}}_4^{\circ }& {\mathit{c}}_5^{\circ }& {\mathit{c}}_6^{\circ }\end{array}\right],$$

$${\mathit{A}}_{\mathsf{\Theta }}^{•}(i,i)=-\mathrm{diag}\left[\begin{array}{cccccc}3.4138& 4.3329& 6.6092& 2.0909& 1& 1\end{array}\right],$$

$${\mathit{A}}_{\mathsf{\Theta }}^{•}(i,i+1)=\mathrm{diag}\left[\begin{array}{cccccc}2.2301& 2.0705& 2.9109& 0& 0& 0\end{array}\right],$$

$${\mathit{A}}_{\mathsf{\Theta }}^{•}(i,i+2)=\mathrm{diag}\left[\begin{array}{cccccc}4.7622& 0.6818& 0& 0& 1& 1\end{array}\right],$$

$${\mathit{A}}_{\mathsf{\Theta }}^{•}(i,i+3)=\mathrm{diag}\left[\begin{array}{cccccc}2.7028& 0& 0& 1.1544& 0& 0\end{array}\right],$$

$${\mathit{A}}_{\mathsf{\Theta }}^{•}(i,i+4)=\mathrm{diag}\left[\begin{array}{cccccc}0& 0& 2.0868& 2.2432& 0& 0\end{array}\right].$$

$${\mathit{A}}_{\mathsf{\Theta }}^{•}(i,i+5)=\mathrm{diag}\left[\begin{array}{cccccc}0& 1.8998& 4.2126& 1.3732& 0& 0\end{array}\right].$$

and the standard parameters that condition the desired specification on the design are set as

$$\mathit{L}=\left[\begin{array}{cc}{0}^{\mathrm{T}}& 1\\ {\mathit{I}}_5& 0\end{array}\right],\mathit{S}=\mathit{L}\otimes {\mathit{I}}_4,{\mathit{J}}^{\mathrm{T}}=\left[\begin{array}{cccccc}{\mathit{I}}_4& {\mathit{I}}_4& {\mathit{I}}_4& {\mathit{I}}_4& {\mathit{I}}_4& {\mathit{I}}_4\end{array}\right].$$

The solution of (55)–(59), obtained using the SeDuMi package [24] in the Matlab environment, is represented by the set of matrix variables

$$\mathit{P}=\mathrm{diag}\left[\begin{array}{cccccc}0.1222& 0.0621& 0.0636& 0.0316& 0.9642& 0.9005\end{array}\right],$$

$$\mathit{R}=\left[\begin{array}{cccc}4.6495& 2.2279& -1.4873& -1.5368\\ 4.1173& 5.5149& -1.8031& -1.8720\\ 0& 0& 0& 0\\ 0& 0& 0& 0\end{array}\right],\mathit{H}=\mathrm{diag}\left[\begin{array}{cccc}0.1222& 0.0621& 0.9642& 0.9005\end{array}\right],$$

which imply the PID control law parameters

$${\mathit{K}}^{\circ }=\left[\begin{array}{cc}{\mathit{K}}_P& -{\mathit{K}}_I\\ 0& 0\end{array}\right]=\left[\begin{array}{cccc}38.0328& 35.8501& -1.5425& -1.7065\\ 33.6799& 88.7411& -1.8700& -2.0788\\ 0& 0& 0& 0\\ 0& 0& 0& 0\end{array}\right],$$

$${\mathit{K}}_P=\left[\begin{array}{cc}38.0328& 35.8501\\ 33.6799& 88.7411\end{array}\right],{\mathit{K}}_I=\left[\begin{array}{cc}1.5425& 1.7065\\ 1.8700& 2.0788\end{array}\right].$$

PID control with these parameters and with the above-defined ${\mathit{K}}_P$ stabilizes the closed-loop system. It is translated into the closed-loop system matrix ${\mathit{A}}_c^{•}$ of the structure

$${\mathit{A}}_c^{•}=\left[\begin{array}{cccccc}\hfill -4.9887& 0.7455& 4.7622\hfill & 2.7028& 0.0639& 0.0707\\ \hfill 1.2093& -6.1524& 2.0705\hfill & 0.6818& 0.0383& 0.0426\\ \hfill 0.5802& 0.9816& -6.6092\hfill & 2.9109& 0.0771& 0.0856\\ \hfill 0.1394& 0.3073& 1.3732\hfill & -2.0909& 0.0498& 0.0553\\ \hfill 1& 0& 0\hfill & 0& -1& 0\\ \hfill 0& 1& 0\hfill & 0& 0& -1\end{array}\right].$$

which is purely Metzler and Hurwitz, with the stable eigenvalue spectrum

$$\rho \left({\mathit{A}}_c^{•}\right)=\left\{\begin{array}{cccccc}-0.3789& -1.0000& -1.2251& -4.2814& -6.8640& -8.0918\end{array}\right\}.$$

Although the matrix ${\mathit{A}}_c^{•}$ is not diagonally dominant, its structure and eigenvalues guarantee that with these PID controller parameters and given non-negative system parameters B, C, a positive closed-loop system performance is achieved. The purpose of the example is primarily to illustrate the desired design procedure.

In

Table 1. The designed parameters dependency on the value of tuning parameter $\epsilon$.

	$\mathit{ϵ}=0.001$	$\mathit{ϵ}=0.01$	$\mathit{ϵ}=0.1$
${\mathit{K}}_D$	$\left[\begin{array}{cc}0.0244& 0\\ 0& 0.0492\end{array}\right]$	$\left[\begin{array}{cc}0.2415& 0\\ 0& 0.4877\end{array}\right]$	$\left[\begin{array}{cc}2.2173& 0\\ 0& 4.4783\end{array}\right]$
${\mathit{K}}_P$	$\left[\begin{array}{cc}38.1003& 36.0151\\ 33.3901& 90.5284\end{array}\right]$	$\left[\begin{array}{cc}38.0328& 35.8501\\ 33.6799& 88.7411\end{array}\right]$	$\left[\begin{array}{cc}37.6024& 34.0588\\ 35.8858& 73.3742\end{array}\right]$
${\mathit{K}}_I$	$\left[\begin{array}{cc}1.5693& 1.7415\\ 1.9075& 2.1252\end{array}\right]$	$\left[\begin{array}{cc}1.5425& 1.7065\\ 1.8700& 2.0788\end{array}\right]$	$\left[\begin{array}{cc}1.3098& 1.4182\\ 1.5143& 1.6492\end{array}\right]$
$\rho \left({\mathit{A}}_c^{•}\right)$	$\left\{\begin{array}{cc}-0.3876,& -1.0000\\ -1.2304,& -4.2752\\ -6.8264,& -8.1135\end{array}\right\}$	$\left\{\begin{array}{cc}-0.3789,& -1.0000\\ -1.2251,& -4.2814\\ -6.8640,& -8.0918\end{array}\right\}$	$\left\{\begin{array}{cc}-0.2826,& -1.0000\\ -1.1755,& -4.3410\\ -7.2581,& -7.8932\end{array}\right\}$

, the main designed parameters presented are dependent on the value of the tuning parameter $\epsilon$.

It is obvious that by increasing value of $\epsilon$, the dominant eigenvalue of the closed-loop structure in the complex plane of the eigenvalues is closer to the imaginary axis, i.e., the dominant time constant is larger. Given all the parameters of the PID controller, the value of tuning parameter $\epsilon =0.01$ is a very good compromise. Moreover, it can be verified that for $\epsilon >0.3$, the closed-loop system is unstable. Using the function $[Y,T,X]=step\left(SYS\right)$ of the Matlab environment, one can easily verify that the closed-loop system state trajectory is aperiodic.

Purely real negative eigenvalues are conditional on the use of positive systems because they guarantee aperiodic positive trajectories of state variables for a non-negative initial state of the system. However, they do not guarantee an overshoot during their evolutions, which sometimes needs to be suppressed. Unfortunately, standard methods for tuning PID controller parameters [25,26] for these structures cannot be used for overshoot suppressing. Methods based on the principle of the D-stability circle region [27] come into consideration, but due to the bilinear structure of the synthesis conditions, it is not possible to guarantee an optimal overshoot suppressing also by using this approach. This sub-area of the synthesis problems will therefore be preferred in authors’ future research.

5. Discussion

The matching conditions restrict the structure of the input system matrix and measured system output variables into a hard constraint on the Metzler non-negative system parameters and limit the PID MIMO implementation. If the matching conditions are satisfied, the D-parameter of the defined PID control law can be tuned to support the Metzler structure of the system matrix, as it can be seen comparing ${\mathit{A}}^{•}$ and A. Considering such support and the feasibility of LMIs design, a closed-loop with a Metzler structure of the dynamics matrix can be expected using the PID controller structure.

Such synthesis task has in general many degrees of freedom in defining structures of the ${\mathit{C}}_p$ matrix. For a given Metzler MIMO system, $\mathit{C}={\mathit{C}}_p$ means the necessary minimal number of measured state variables, which is not changed if any of the necessary measured state variables for the implementation of the D-component part of PID are not used in the construction of the I-component part. The structure of the matrix $\mathit{C}={\mathit{C}}_p$ was chosen to demonstrate solutions with equal measured state impact also on the I-component part. For the above-considered structure of B, an additional measured state variable can be included only in the I-component part but the matching conditions imply that a strictly Metzler ${\mathit{A}}_c^{•}$ cannot be obtained.

To the best of the authors’ knowledge, no comparable results are available for the design of PID control of MIMO Metzler linear systems. In the authors’ opinion, the proposed method is one that gives constraint limits on conditions for a class of switched positive systems. Exposing the principle details, the approach can be adapted for studying the PI control of strictly Metzler MIMO linear systems, where similar results can be expected.

6. Conclusions

This paper completes a design method for the synthesis of PID control for Metzler continuous-time MIMO linear systems. The closed-loop purely Metzler system matrix structure is proposed when exploiting a tuned diagonal D-part gain matrix and the square positive D-part and I-part gain matrices. The newly formulated exposition of the problem treatments of the existing matching conditions were provided by a measurement assignment through the output matrix structure to find an LMI representation of the design conditions. Maintaining parametric system constraints by the set of LMIs, the design conditions were completed by Lyapunov matrix inequality, guaranteeing closed-loop asymptotic stability within a feasible solution.

Since the analysis is linear, evidently, one can see the dependence of the resulting PID gains on Metzler parameters of the system. The proposed approach lends itself to algorithm formalization through LMIs. The theory yields results that have otherwise not been derived for PID control of a given class of Metzler systems. The development of an approach for Metzler systems with an extended set of parametric constraints is a topic of future research.