A Randomized Algorithm for Optimal PID Controllers

Abstract

A randomized algorithm is suggested for the syntheses of optimal PID controllers for MIMO coupled systems, where the optimality is with respect to the $H_{\infty }$-norm, the $H_2$-norm and the LQR functional, with possible system-performance specifications defined by regional pole-placement. Other notions of optimality (e.g., mixed $H_2/H_{\infty }$ design, controller norm or controller sparsity) can be handled similarly with the suggested algorithm. The suggested method is direct and thus can be applied to continuous-time systems as well as to discrete-time systems with the obvious minor changes. The presented algorithm is a randomized algorithm, which has a proof of convergence (in probability) to a global optimum.

1. Introduction

Proportional Integral Derivative (PID) controllers are the most widely-used controllers in industry despite of many new results in control theory that have been achieved in recent years. This is because in many cases PID controllers achieve close to optimal performance, because it is hard to compete with their investment-to-profit ratio and because of their well-understood influence on the closed-loop system behavior, as well as the availability of many tuning procedures for various performance objectives, at least for SISO systems or MIMO systems that can be well decoupled. Often, the performance achieved by PID controllers can be significantly improved by using more effective design techniques. Therefore the field of PID controllers is an on-going active field of research. A review of most recent techniques for tuning and designing PID-based control structures, together with methods for assessing their performance is presented in [1]. In [2] optimal decentralized MIMO PID controller design for a benchmark MIMO distillation column models is carried out using a Modified Firefly Algorithm (MFA), where the optimality is with respect to an objective function which is a weighted combination of the $L_1$-norm and the $L_2$-norm of the tracking error, with respect to a reference signal.

Tuning of MIMO coupled systems PID controllers is much more challenging than tuning PID controllers for SISO system or decoupled MIMO systems, because of the number of parameters and because they should be tuned all at once, since otherwise, tuning each channel separately might end up in a very poor performance of the whole system or even destabilize the whole system. In [3] the need for effective methods for the design of PID optimal controllers for MIMO coupled systems that cannot be well approximated by second-order systems, was raised. In [4] the parametrization of all the solutions of the standard $H_{\infty }$ control problem is used to obtain Bilinear Matrix Inequalities (BMI’s) constraints on the PID gains for each frequency. Next, for each frequency, Linear Matrix Inequalities (LMI’s) restricted constraints for a stabilizing PID gain are derived from the BMI constraints. Finally, the set of LMI’s for a chosen finite set of frequencies are solved by the Iterated Linear Matrix Inequalities (ILMI) method. The sequence is shown to converge to a local optimum. In [5] a method for MIMO coupled systems PID controller design for stable plants (given by transfer function at an appropriate set of frequencies) is represented. The method relies on solving a sequence of SDP’s, but it cannot guarantee finding a globally optimal design.

In this article, we suggest a randomized algorithm for PID optimal controllers. The algorithm is based on the known reduction of [3] to the Stabilizing Static-Output-Feedback (SOF) problem, which is then solved by the Ray-Shooting Method, where the method is used both for the search for feasible starting SOF and for the optimal SOF (see [6,7,8] for some applications of the method).

A closely related approach to ours is given in [9] where the problem is first reduced to the SOF problem, by the reduction of [3] and by using a descriptor system representation of the augmented system (in order to bypass the feasibility condition that $I-{\widehat{K}}_3{C}_2{B}_2$ should be invertible—see the following). Next, the search problem of a starting stabilizing SOF controller is solved via the ILMI approach. Finally, a Bounded-Real Lemma version for descriptor systems is used in order to get a set of Quadratic Matrix Inequalities (QMI’s) that their solution guaranties a $\gamma$ suboptimal $H_{\infty }$ gain. The last can be solved by the ILMI method (without a proof of convergence). The method of ILMI was discussed in many articles (see [3,10,11,12]). Unfortunately, the ILMI method needs a good starting point and has no proof of convergence in general.

Another closely related approach is given in [13] where the 2-Degree of Freedom PID (2-DOF-PID) problem is reduced to the SOF problem, which in order to incorporate the optimization problem of $\gamma$ suboptimal $H_{\infty }$ gain, is reformulated as a set of BMI’s, where the last is solved via standard BMI solvers.

Many problems can be reduced to the SOF problem, imposing structural or other constraints on the controllers. These include the reduction of the minimal-degree dynamic-feedback problem and robust or decentralized stability via static-feedback (see [14,15], resp.). The design problem of remote digital PID controllers formulated as robust SOF problem for discrete-time Networked Control Systems (NCS), subject to time-delays and randomly missing-data is considered in [16]. A formulation of the reduced-order $H_{\infty }$ filter problem as constrained SOF problem, is considered in [17]. Unfortunately, the SOF problem with bounded entries or other structural constraints is known to be NP-hard (see [14,18]). Therefore, minimal-gain SOF with bounded entries is obviously NP-hard problem. Exact Pole-placement and simultaneous and robust stabilization via SOF are also NP-hard (see [14,19], resp.). Thus, practically, one should expect that only approximation or randomized algorithms will be able to cope with these problems.

The article is organized as follows:

In Section 2 we define the feasibility problem of PID controllers and the specific reduction to the SOF problem that will be used. In Section 3 we define the PID controller optimization problems where the optimality is with respect to the $H_{\infty }$-norm, the $H_2$-norm and the LQR functional, with possible regional pole-placement. In Section 4 we describe the suggested randomized algorithm and finally, in Section 5 we give a real-life example for each one of the above-mentioned optimization problems.

The notions are standard. By $\Re \left(z\right)$ and $\Im \left(z\right)$ we denote the real and imaginary parts of $z\in \mathbb{C}$, resp., where by ${\mathbb{C}}_{-}$ we denote the left half-plane. For a matrix Z we denote by $Z^T$ its transpose and by $\overline{Z}$ we denote the conjugate matrix. By $Z^{*}$ we denote conjugate-transpose of Z. By $Z^{+}$ we denote the Moore-Penrose pseudo-inverse of Z, where by $L_Z$ and $R_Z$ we denote the left and right orthogonal projections $I-Z^{+}Z$ and $I-Z{Z}^{+}$, resp.. For a square matrix Z, we denote by $\sigma \left(Z\right)$ the spectrum of Z and by $trace\left(Z\right)$ we denote the trace of Z. For matrices $Z,W$ we denote by $Z\otimes W$ the Kronecker product. By $∥Z∥$ we denote the spectral norm of Z, i.e., the largest singular value of Z, where by ${∥Z∥}_F$ we denote its Frobenious norm, i.e., ${\left(trace\left(Z^{*}Z\right)\right)}^{\frac{1}{2}}$. For matrices $Z,W$ with the same size, we denote by ${〈Z,W〉}_F=trace\left(W^{*}Z\right)$ the Frobenious inner product. For a vector v we denote by $∥v∥$ its Euclidean norm ${\left({\sum }_j{\left|v_j\right|}^2\right)}^{\frac{1}{2}}$. By $vec\left(Z\right)$ we denote the vector obtained from the matrix Z by chaining all its columns to a single column. By $mat$ we denote the inverse function of $vec$. For a transfer function $T\left(s\right)$ analytic in the open right half-plane, we denote by ${∥T∥}_{H_{\infty }}$ its $H_{\infty }$ norm where ${∥T∥}_{H_{\infty }}=ess-su{p}_{\omega \in \mathbb{R}}∥T\left(j\omega \right)∥$ and by ${∥T∥}_{H_2}$ we denote its $H_2$ norm where ${∥T∥}_{H_2}={\left(\frac{1}{2\pi }{\int }_{-\infty }^{+\infty }trace\left(T{\left(j\omega \right)}^{*}T\left(j\omega \right)\right)d\omega \right)}^{\frac{1}{2}}$. For a set X in a topological space, we denote by $\overline{X}$ its closure and by $int\left(X\right)$ its interior. For two sets $X,Y$ we denote by $CH\left(X,Y\right)$ the convex-hull of X and Y, i.e., the minimal closed convex set that contains $X\cup Y$.

2. The Problem of Stabilization via PID Controllers

Let a system be given by:

$$\left\{\begin{array}{c}\begin{array}{cc}\hfill \dot{x}\left(t\right)& =Ax\left(t\right)+B_{1}w\left(t\right)+B_{2}u\left(t\right)\hfill \\ \hfill z\left(t\right)& =C_{1}x\left(t\right)+D_{1,1}w\left(t\right)+D_{1,2}u\left(t\right)\hfill \\ \hfill y\left(t\right)& =C_{2}x\left(t\right)+D_{2,1}w\left(t\right),\hfill \end{array}\hfill \end{array}\right.$$

where $A\in {\mathbb{R}}^{n_x\times n_x},B_1\in {\mathbb{R}}^{n_x\times n_w},B_2\in {\mathbb{R}}^{n_x\times n_u},C_1\in {\mathbb{R}}^{n_z\times n_x},D_{1,1}\in {\mathbb{R}}^{n_z\times n_w},D_{1,2}\in {\mathbb{R}}^{n_z\times n_u},C_2\in {\mathbb{R}}^{n_y\times n_x},D_{2,1}\in {\mathbb{R}}^{n_y\times n_w}$, where x is the state, z is the regulated output, y is the measurement, u is the control input and w is the noise. Assume that $D_{2,1}=0$ and let

$$u\left(t\right)=r\left(t\right)-K_{P}y\left(t\right)-K_I{\int }_0^{t}y\left(\tau \right)d\tau -K_D\dot{y}\left(t\right),$$

where $r\left(t\right)$ is the new reference input. Let $\widehat{x}\left(t\right)=\left[\begin{array}{c}x\left(t\right)\hfill \\ {\int }_0^{t}y\left(\tau \right)d\tau \hfill \end{array}\right]$ be the augmented state, then:

$$\begin{array}{cc}\hfill u\left(t\right)& =r\left(t\right)-K_{P}y\left(t\right)-K_I{\int }_0^{t}y\left(\tau \right)d\tau -K_D\dot{y}\left(t\right)\hfill \\ & =r-K_P\left[\begin{array}{cc}{C}_2& 0\end{array}\right]\widehat{x}\left(t\right)-K_I\left[\begin{array}{cc}0& I\end{array}\right]\widehat{x}\left(t\right)+\hfill \\ & -K_D{C}_2\dot{x}\left(t\right)\hfill \\ & =r-K_P\left[\begin{array}{cc}{C}_2& 0\end{array}\right]\widehat{x}\left(t\right)-K_I\left[\begin{array}{cc}0& I\end{array}\right]\widehat{x}\left(t\right)+\hfill \\ & -K_D{C}_2\left(Ax\left(t\right)+B_{1}w\left(t\right)+B_{2}u\left(t\right)\right),\hfill \end{array}$$

implying that:

$$\begin{array}{cc}\hfill u\left(t\right)& =M^{-1}r\left(t\right)-M^{-1}{K}_P\left[\begin{array}{cc}{C}_2& 0\end{array}\right]\widehat{x}\left(t\right)+\hfill \\ & -M^{-1}{K}_I\left[\begin{array}{cc}0& I\end{array}\right]\widehat{x}\left(t\right)+\hfill \\ & -M^{-1}{K}_D\left[\begin{array}{cc}{C}_{2}A& 0\end{array}\right]\widehat{x}\left(t\right)-M^{-1}{K}_D{C}_2{B}_{1}w,\hfill \end{array}$$

where $M=I+K_D{C}_2{B}_2$ is assumed to be invertible. Now:

$$\begin{array}{cc}\hfill \dot{\widehat{x}}\left(t\right)& =\left[\begin{array}{c}\dot{x}\left(t\right)\\ y\left(t\right)\end{array}\right]\hfill \\ & =\left[\begin{array}{c}Ax\left(t\right)+B_{1}w\left(t\right)+B_{2}u\left(t\right)\\ C_{2}x\left(t\right)\end{array}\right]\hfill \\ & =\left[\begin{array}{cc}A& 0\\ C_2& 0\end{array}\right]\widehat{x}\left(t\right)+\hfill \\ & +\left[\begin{array}{c}{B}_1-B_2{M}^{-1}{K}_D{C}_2{B}_1\\ 0\end{array}\right]\widehat{w}\left(t\right)+\left[\begin{array}{c}{B}_2\\ 0\end{array}\right]\widehat{u}\left(t\right),\hfill \end{array}$$

where:

$$\left\{\begin{array}{c}\widehat{w}\left(t\right)=w\left(t\right)\hfill \\ \widehat{u}\left(t\right)=\widehat{r}\left(t\right)-\widehat{K}\widehat{y}\left(t\right)\hfill \\ \widehat{r}\left(t\right)=M^{-1}r\left(t\right)\hfill \\ \widehat{y}\left(t\right)=\left[\begin{array}{c}{C}_{2}x\left(t\right)\\ {\int }_0^{t}y\left(\tau \right)d\tau \\ C_{2}Ax\left(t\right)\end{array}\right]\hfill \\ \widehat{K}=M^{-1}\left[\begin{array}{ccc}{K}_P& K_I& K_D\end{array}\right].\hfill \end{array}\right.$$

Let

$$\left\{\begin{array}{c}\widehat{A}=\left[\begin{array}{cc}A& 0\\ C_2& 0\end{array}\right],{\widehat{B}}_2=\left[\begin{array}{c}{B}_2\\ 0\end{array}\right]\hfill \\ {\widehat{B}}_1=\left[\begin{array}{c}{B}_1-B_2{M}^{-1}{K}_D{C}_2{B}_1\\ 0\end{array}\right]\hfill \\ {\widehat{C}}_1=\left[\begin{array}{cc}{C}_1& 0\end{array}\right],{\widehat{C}}_2=\left[\begin{array}{cc}{C}_2& 0\\ 0& I\\ C_{2}A& 0\end{array}\right]\hfill \\ {\widehat{D}}_{1,1}=D_{1,1}-D_{1,2}{M}^{-1}{K}_D{C}_2{B}_1,{\widehat{D}}_{1,2}=D_{1,2}.\hfill \end{array}\right.$$

Let $\widehat{z}\left(t\right)=z\left(t\right)$ then, the augmented system is given by:

$$\left\{\begin{array}{c}\dot{\widehat{x}}\left(t\right)=\left(\widehat{A}-{\widehat{B}}_2\widehat{K}{\widehat{C}}_2\right)\widehat{x}\left(t\right)+{\widehat{B}}_1\widehat{w}\left(t\right)+{\widehat{B}}_2\widehat{r}\left(t\right)\hfill \\ \widehat{z}\left(t\right)=\left({\widehat{C}}_1-{\widehat{D}}_{1,2}\widehat{K}{\widehat{C}}_2\right)\widehat{x}\left(t\right)+{\widehat{D}}_{1,1}\widehat{w}\left(t\right)+{\widehat{D}}_{1,2}\widehat{r}\left(t\right)\hfill \\ \widehat{y}\left(t\right)={\widehat{C}}_2\widehat{x}\left(t\right).\hfill \end{array}\right.$$

(1)

Once that a stabilizing SOF matrix $\widehat{K}=\left[\begin{array}{ccc}{\widehat{K}}_1\hfill & {\widehat{K}}_2\hfill & {\widehat{K}}_3\hfill \end{array}\right]$ for the augmented system (1) has been found, we have $K_D=M{\widehat{K}}_3$. Thus, $M=I+M{\widehat{K}}_3{C}_2{B}_2$ from which we conclude that:

$$M={\left(I-{\widehat{K}}_3{C}_2{B}_2\right)}^{-1}.$$

Finally, we have:

$$\left\{\begin{array}{c}{K}_P=M{\widehat{K}}_1={\left(I-{\widehat{K}}_3{C}_2{B}_2\right)}^{-1}{\widehat{K}}_1\hfill \\ K_I=M{\widehat{K}}_2={\left(I-{\widehat{K}}_3{C}_2{B}_2\right)}^{-1}{\widehat{K}}_2\hfill \\ K_D=M{\widehat{K}}_3={\left(I-{\widehat{K}}_3{C}_2{B}_2\right)}^{-1}{\widehat{K}}_3.\hfill \end{array}\right.$$

(2)

We conclude that the stabilization problem via PID controller is reducible to the stabilization via SOF problem, with the constraint that $I-{\widehat{K}}_3{C}_2{B}_2$ is invertible. In the following we will assume that the couples $\left(\widehat{A},{\widehat{B}}_2\right)$ and $\left({\widehat{A}}^T,{\widehat{C}}_2^T\right)$ are controllable, in order to make the task of finding a feasible starting SOF $\widehat{K}$ easier.

Remark 1.

Please note that although the method relies on a specific known state-space model of the plant, it does not require an exact model, since the modeling errors can be included in the dynamic uncertainty and they can be taken into account by selecting a nominal model with a suitable weighting function (see [20]). In any case, after designing a PID controller for the nominal model, its performance should be checked on the non-linear plant model. Note also that the requirement that $I-{\widehat{K}}_3{C}_2{B}_2$ should be invertible is not a problem since that if a stabilizing SOF $\widehat{K}$ exists, there exist an open neighborhood of it where all the members are stabilizing SOF’s, where on the other hand, the set of all matrices that does not satisfy the invertibility condition has measure 0.

Remark 2.

Since the problem of SOF is (polynomial-time) reducible to the PID controller problem (by restriction to $K_I=0,K_D=0$), therefore the PID controller problem is as hard as the SOF problem. Since the SOF problem with structural constraints (on $K_P$), or the problem of exact pole-placement via SOF are NP-hard, it follows that these problems with PID controllers are NP-hard. We therefore conclude that unless $P=NP$, no deterministic polynomial-time algorithm can solve these problems and we therefore are led to use randomization with the hope that randomized algorithms will be able to cope more efficiently with these problems.

3. Statement of the Problem

The $H_{\infty }$ and $H_2$ norm minimization relates to the transfer function $T_{\widehat{w},\widehat{z}}\left(s\right)$ from the noise $\widehat{w}$ to the output $\widehat{z}$ (when the reference signal $\widehat{r}$ is zero and the system is derived only by the noise). The transfer function is given by:

$$T_{\widehat{w},\widehat{z}}\left(s\right)={\widehat{D}}_{1,1}+\left({\widehat{C}}_1-{\widehat{D}}_{1,2}\widehat{K}{\widehat{C}}_2\right){\left(sI-\widehat{A}+{\widehat{B}}_2\widehat{K}{\widehat{C}}_2\right)}^{-1}{\widehat{B}}_1.$$

(3)

In order that the $H_2$-norm of $T_{\widehat{w},\widehat{z}}\left(s\right)$ would be finite, we need that ${\widehat{D}}_{1,1}=0$. We therefore have the following lemma (given without a proof):

Lemma 1.

We have:

$${\widehat{D}}_{1,1}=D_{1,1}-D_{1,2}{M}^{-1}{K}_D{C}_2{B}_1=D_{1,1}-D_{1,2}{\widehat{K}}_3{C}_2{B}_1=0,$$

(4)

if and only if:

$$\left({\left(C_2{B}_1\right)}^T\otimes D_{1,2}\right)·{\left({\left(C_2{B}_1\right)}^T\otimes D_{1,2}\right)}^{+}vec\left(D_{1,1}\right)=vec\left(D_{1,1}\right).$$

(5)

In this case, the set of all matrices ${\widehat{K}}_3$ satisfying (4) is given by:

$${\widehat{K}}_3=mat\left({\left({\left(C_2{B}_1\right)}^T\otimes D_{1,2}\right)}^{+}vec\left(D_{1,1}\right)+L_{{\left(C_2{B}_1\right)}^T\otimes D_{1,2}}vec\left(Z\right)\right),$$

(6)

where Z is arbitrary matrix with the same size as ${\widehat{K}}_3$.

For the LQR minimization, we assume that $\widehat{w}\left(t\right)=0$ and $\widehat{r}\left(t\right)=0$, i.e., the system is derived by the initial state ${\widehat{x}}_0:=\widehat{x}\left(0\right)$ (initial step input). Assume that a SOF $\widehat{K}$ was found, such that ${\widehat{A}}_{c\ell }\left(\widehat{K}\right):=\widehat{A}-{\widehat{B}}_2\widehat{K}{\widehat{C}}_2$ is stable. Let the LQR functional be defined by:

$$J\left({\widehat{x}}_0\right):={\int }_0^{\infty }\left(\widehat{x}{\left(t\right)}^{T}Q\widehat{x}\left(t\right)+\widehat{z}{\left(t\right)}^{T}R\widehat{z}\left(t\right)\right)dt,$$

(7)

where $Q>0$ and $R\ge 0$. Let $\widehat{C}\left(\widehat{K}\right):={\widehat{C}}_1-{\widehat{D}}_{1,2}\widehat{K}{\widehat{C}}_2$. Then, $\widehat{z}\left(t\right)=\widehat{C}\left(\widehat{K}\right)\widehat{x}\left(t\right)$ and substitution of the last into (7) yields

$$J\left({\widehat{x}}_0,\widehat{K}\right)={\int }_0^{\infty }\widehat{x}{\left(t\right)}^T\left(Q+\widehat{C}{\left(\widehat{K}\right)}^{T}R\widehat{C}\left(\widehat{K}\right)\right)\widehat{x}\left(t\right)dt={\widehat{x}}_0^{T}P\left(\widehat{K}\right){\widehat{x}}_0,$$

(8)

where $P\left(\widehat{K}\right)>0$ is the unique solution to the Lyapunov equation

$${\widehat{A}}_{c\ell }{\left(\widehat{K}\right)}^{T}P+P{\widehat{A}}_{c\ell }\left(\widehat{K}\right)=-\left(Q+\widehat{C}{\left(\widehat{K}\right)}^{T}R\widehat{C}\left(\widehat{K}\right)\right)$$

(9)

given explicitly by

$$P\left(\widehat{K}\right)=-mat\left({\left(I\otimes {\widehat{A}}_{c\ell }{\left(\widehat{K}\right)}^T+{\widehat{A}}_{c\ell }{\left(\widehat{K}\right)}^T\otimes I\right)}^{-1}·vec\left(Q+\widehat{C}{\left(\widehat{K}\right)}^{T}R\widehat{C}\left(\widehat{K}\right)\right)\right).$$

(10)

Thus, we look for SOF $\widehat{K}$ that minimizes the functional $J\left(x_0,\widehat{K}\right)={\widehat{x}}_0^{T}P\left(\widehat{K}\right){\widehat{x}}_0$. When ${\widehat{x}}_0$ is unknown, we seek a SOF $\widehat{K}$ for which

$${\sigma }_{max}\left(\widehat{K}\right):=max\left(\sigma \left(P\left(\widehat{K}\right)\right)\right)$$

(11)

is minimal. In this case, we get a robust LQR via SOF, in the sense that it minimizes $J\left({\widehat{x}}_0,\widehat{K}\right)$ for the worst possible (unknown) ${\widehat{x}}_0$. Please note that

$$\begin{array}{cc}\hfill {\widehat{x}}_0^{T}P\left(\widehat{K}\right){\widehat{x}}_0& ={∥P{\left(\widehat{K}\right)}^{\frac{1}{2}}{\widehat{x}}_0∥}^2\le \hfill \\ & \le {∥P{\left(\widehat{K}\right)}^{\frac{1}{2}}∥}^2{∥{\widehat{x}}_0∥}^2=\hfill \\ & =∥P\left(\widehat{K}\right)∥{∥{\widehat{x}}_0∥}^2={\sigma }_{max}\left(\widehat{K}\right){∥{\widehat{x}}_0∥}^2,\hfill \end{array}$$

and that there exists ${\widehat{x}}_0\ne 0$ for which equality holds. Therefore $\frac{J\left({\widehat{x}}_0,\widehat{K}\right)}{{∥{\widehat{x}}_0∥}^2}\le {\sigma }_{max}\left(\widehat{K}\right)$, with equality in the worst case.

Let $\Omega$ denote a subset of ${\mathbb{C}}_{-}$, which is symmetric with respect to the x-axis, containing a positive length segment of the real axis, with some neighborhood of it, and such that $\Omega =\overline{int\left(\Omega \right)}$ (i.e., $\Omega$ is the closure of its interior). Let ${\mathcal{S}}^{q\times r}$ denote the set of all matrices $\widehat{K}\in {\mathbb{R}}^{q\times r}$, such that ${\widehat{A}}_{c\ell }\left(\widehat{K}\right)$ is stable, i.e., $\sigma \left({\widehat{A}}_{c\ell }\left(\widehat{K}\right)\right)\subset {\mathbb{C}}_{-}$. By ${\mathcal{S}}_{\Omega }^{q\times r}$, we denote the set of all matrices $\widehat{K}\in {\mathbb{R}}^{q\times r}$, such that $\sigma \left({\widehat{A}}_{c\ell }\left(\widehat{K}\right)\right)\subset \Omega$. In this case, we say that ${\widehat{A}}_{c\ell }\left(\widehat{K}\right)$ is $\Omega$-stable (or, the augmented system (1) is $\Omega$-stable). In the sequel, we will occasionally write ${\mathcal{S}}_{\Omega }$ instead of ${\mathcal{S}}_{\Omega }^{q\times r}$, when it is clear what the size of the related matrices is. Please note that ${\mathcal{S}}_{\Omega }$ is a closed set, since $\Omega$ is closed and that ${\mathcal{S}}_{\Omega }$ has nonempty interior (if it is nonempty) since $\Omega$ has nonempty interior.

Under the assumptions that $D_{2,1}=0$ and that $\left(\widehat{A},{\widehat{B}}_2\right),\left({\widehat{A}}^T,{\widehat{C}}_2^T\right)$ are controllable (and the additional assumption (5) when dealing with the $H_2$-norm minimization), for the PID controller problem with minimal $H_{\infty }$-norm, we need to solve:

$$\begin{array}{c}min{f}_{H_{\infty }}\left(\widehat{K}\right):={∥T_{\widehat{w},\widehat{z}}∥}_{H_{\infty }}\hfill \\ \mathrm{such}\mathrm{that}:\hfill \\ \left\{\begin{array}{c}\widehat{K}\in {\mathcal{S}}_{\Omega }\hfill \\ I-{\widehat{K}}_3{C}_2{B}_2\mathrm{is}\mathrm{invertible}\hfill \end{array}\right.\hfill \end{array}$$

(12)

where $T_{\widehat{w},\widehat{z}}$ is given by (3). For the PID controller problem with minimal $H_2$-norm, we need to solve:

$$\begin{array}{c}min{f}_{H_2}\left(\widehat{K}\right):={∥T_{\widehat{w},\widehat{z}}∥}_{H_2}\hfill \\ \mathrm{such}\mathrm{that}:\hfill \\ \left\{\begin{array}{c}\widehat{K}\in {\mathcal{S}}_{\Omega }\hfill \\ {\widehat{K}}_3\mathrm{has}\mathrm{the}\mathrm{structure}\left(6\right)\hfill \\ I-{\widehat{K}}_3{C}_2{B}_2\mathrm{is}\mathrm{invertible}\hfill \end{array}\right.\hfill \end{array}$$

(13)

where $T_{\widehat{w},\widehat{z}}$ is given by (3) and, for the PID controller problem with minimal LQR functional-value, we need to solve:

$$\begin{array}{c}min{f}_{LQR}\left(\widehat{K}\right):={\sigma }_{max}\left(\widehat{K}\right)\hfill \\ \mathrm{such}\mathrm{that}:\hfill \\ \left\{\begin{array}{c}\widehat{K}\in {\mathcal{S}}_{\Omega }\hfill \\ I-{\widehat{K}}_3{C}_2{B}_2\mathrm{is}\mathrm{invertible}\hfill \end{array}\right.\hfill \end{array}$$

(14)

where ${\sigma }_{max}\left(\widehat{K}\right)$ is given by (11).

4. The Suggested Algorithm

Assume that ${\widehat{K}}^{\left(0\right)}\in int{\left({\mathcal{S}}_{\Omega }\right)}^{q\times r}$ was found by the Ray-Shooting algorithm ([6]) or by any other method (see [21,22,23]). We assume that $I-{\widehat{K}}_3^{\left(0\right)}{C}_2{B}_2$ is invertible (and in addition, we assume that ${\widehat{K}}_3^{\left(0\right)}$ has the form (6) for some initial matrix $Z^{\left(0\right)}$, when dealing with the $H_2$-norm minimization). Let $h>0$ and let $U^{\left(0\right)}$ be a unit vector w.r.t. the Frobenius norm, i.e., ${∥U^{\left(0\right)}∥}_F=1$. Let $L^{\left(0\right)}={\widehat{K}}^{\left(0\right)}+h·U^{\left(0\right)}$ and let $\mathcal{L}$ be the hyperplane defined by $L^{\left(0\right)}+V$, where ${〈V,U^{\left(0\right)}〉}_F=0$. Let $r_{\infty }>0$ and let ${\mathcal{R}}_{\infty }$ denote the set of all $F\in \mathcal{L}$, such that ${∥F-L^{\left(0\right)}∥}_F\le r_{\infty }$. Let ${\mathcal{R}}_{\infty }\left(ϵ\right)={\mathcal{R}}_{\infty }+\overline{\mathbb{B}\left(0,ϵ\right)}$, where $\overline{\mathbb{B}\left(0,ϵ\right)}$ denotes the closed ball centered at 0 with radius $ϵ$ ($0<ϵ\le \frac{1}{2}$), with respect to the Frobenius norm on ${\mathbb{R}}^{q\times r}$. Let ${\mathcal{D}}^{\left(0\right)}=CH\left(K^{\left(0\right)},{\mathcal{R}}_{\infty }\left(ϵ\right)\right)$ denote the convex-hull of the vertex ${\widehat{K}}^{\left(0\right)}$ with the basis ${\mathcal{R}}_{\infty }\left(ϵ\right)$. Let ${\mathcal{S}}_{\Omega }^{\left(0\right)}={\mathcal{S}}_{\Omega }\cap {\mathcal{D}}^{\left(0\right)}$ and please note that ${\mathcal{S}}_{\Omega }^{\left(0\right)}$ is compact (but generally not convex). We wish to minimize a continuous function $f\left(\widehat{K}\right)$ over the compact set ${\mathcal{S}}_{\Omega }\cap \overline{\mathbb{B}\left({\widehat{K}}^{\left(0\right)},h\right)}$. Let ${\widehat{K}}_{*}$ denote a point in ${\mathcal{S}}_{\Omega }\cap \overline{\mathbb{B}\left({\widehat{K}}^{\left(0\right)},h\right)}$ where a global minimum of $f\left(\widehat{K}\right)$ is accepted. Obviously, ${\widehat{K}}_{*}\in {\mathcal{D}}^{\left(0\right)}$, for some direction $U^{\left(0\right)}$, as above. When dealing with the $H_2$-norm minimization, and if $U^{\left(0\right)}=\left[\begin{array}{ccc}{U}_1^{\left(0\right)}& U_2^{\left(0\right)}& U_3^{\left(0\right)}\end{array}\right]$ then, the directions will be constrained to have the form $U_3^{\left(0\right)}=mat\left(L_{{\left(C_2{B}_1\right)}^T\otimes D_{1,2}}vec\left(Z\right)\right)$, in order that $L^{\left(0\right)}$ would have the form (6). We will also restrict V as we have restricted $U^{\left(0\right)}$, in order that any element of $\mathcal{L}$ would have the structure (6).

The suggested Algorithm 1 goes as follows:

We start with a point ${\widehat{K}}^{\left(0\right)}\in int\left({\mathcal{S}}_{\Omega }\right)$, found by the Ray-Shooting algorithm [6]. Assuming that ${\widehat{K}}_{*}\in {\mathcal{D}}^{\left(0\right)}$, the inner-loop ($j=1,\dots ,n$) uses the Ray-Shooting Method in order to find an approximation of the global minimum of the function $f\left(\widehat{K}\right)$ over ${\mathcal{S}}_{\Omega }^{\left(0\right)}$ - the portion of ${\mathcal{S}}_{\Omega }$ bounded in the cone ${\mathcal{D}}^{\left(0\right)}$. The proof of convergence in probability of the inner-loop and its complexity can be found in [6] (see also [7]). In the inner-loop, we choose a search direction by choosing a point F in ${\mathcal{R}}_{\infty }\left(ϵ\right)$ - the base of the cone ${\mathcal{D}}^{\left(0\right)}$. Next, in the most inner-loop ($k=1,\dots ,s$) we scan the ray $\widehat{K}\left(t\right):=\left(1-t\right){\widehat{K}}^{\left(0\right)}+tF$ and record the best controller on it. Repeating this sufficiently many times, we reach ${\widehat{K}}_{*}$ (or an $ϵ$- neighborhood of it) with high probability, under the assumption that ${\widehat{K}}_{*}\in {\mathcal{D}}^{\left(0\right)}$. In [8] it was shown that by taking $h=r_{\infty }$ and $m=⌈e·\sqrt{2\pi \ell }⌉$ iterations in the outer-loop (where $\ell :=qr$), we have ${\widehat{K}}_{*}\in {\mathcal{D}}^{\left(0\right)}$ almost surely. Specifically, when $\ell \ge 12$, it was suggested to take $m=2\ell$ and $\ell \times \ell$ orthogonal matrix $U=\left[\begin{array}{cccc}{u}_1& u_2& \cdots & u_{\ell }\end{array}\right]$, and to take the directions $U_j^{\left(0\right)}=\pm mat\left(u_j\right),j=1,\dots ,\ell$ in the outer-loop.

The outer-loop ($i=1,\dots ,m$) is used instead of executing the Ray-Shooting algorithm again and again, by taking ${\widehat{K}}^{\left(\mathrm{best}\right)}$ as the new vertex of the search cone instead of ${\widehat{K}}^{\left(0\right)}$. The replacement of ${\widehat{K}}^{\left(0\right)}$ can be avoided if $r_{\infty }$ is chosen sufficiently large. The replacement of ${\widehat{K}}^{\left(0\right)}$ by ${\widehat{K}}^{\left(\mathrm{best}\right)}$ can be considered as a heuristic step, which is made instead of running the Ray-Shooting algorithm many times in order to generate "the best starting point", which is relevant only if we actually evaluate $f\left(\widehat{K}\right)$ on each such point and take the point with the best value as the best starting point. Since in any case, we evaluate $f\left(\widehat{K}\right)$ in the main algorithm, we could avoid the repeated execution of the Ray-Shooting algorithm. The outer-loop is similar to what is done in the Hide-And-Seek algorithm (see [24,25]). The convergence in probability of the Hide-And-Seek algorithm can be found in [26].

Algorithm 1 Optimal PID controller randomized algorithm

Require:  An algorithm for deciding Ω-stability and an algorithm for computing $f\left(\widehat{K}\right)$.    Input: ${\widehat{K}}^{\left(0\right)}$ such that ${\widehat{A}}_{c\ell }\left({\widehat{K}}^{\left(0\right)}\right)$ is Ω-stable and $I-{\widehat{K}}_3^{\left(0\right)}{C}_2{B}_2$ is invertible (with the needed structure (6) if $H_2$-norm minimization is sought), $ϵ>0,r_{\infty }>0,h>0$ and integers $m,n,s$.    Output: ${\widehat{K}}_{*}$ that globally minimize $f\left(\widehat{K}\right)$ (or within ϵ distance from a global minimum of $f\left(\widehat{K}\right)$) in h-radius closed-ball centered at ${\widehat{K}}^{\left(0\right)}$.  1. ${\widehat{K}}^{\left(best\right)}\leftarrow {\widehat{K}}^{\left(0\right)}$  2. $f^{\left(best\right)}\leftarrow f\left({\widehat{K}}^{\left(0\right)}\right)$  3. $i\leftarrow 0$  4. $k\leftarrow 1$  5. while$i\le m$do  6.  $i\leftarrow i+1$  7.  choose $U^{\left(0\right)}$ such that ${∥U^{\left(0\right)}∥}_F=1$, uniformly at random  8.  set $L^{\left(0\right)}\leftarrow {\widehat{K}}^{\left(0\right)}+h·U^{\left(0\right)}$  9.  for $j=1$ to n do 10.   choose $F\in {\mathcal{R}}_{\infty }\left(ϵ\right)$ uniformly at random 11.   for $k=1$ to s do 12.    $t\leftarrow \frac{k}{s}$ 13.    $\widehat{K}\left(t\right)\leftarrow \left(1-t\right){\widehat{K}}^{\left(0\right)}+tF$ 14.    if $\widehat{K}\left(t\right)\in {\mathcal{S}}_{\Omega }$ then 15.     if $f\left(\widehat{K}\left(t\right)\right)<f^{\left(best\right)}$ then 16.      ${\widehat{K}}^{\left(best\right)}\leftarrow \widehat{K}\left(t\right)$ 17.      $f^{\left(best\right)}\leftarrow f\left(\widehat{K}\left(t\right)\right)$ 18.     end if 19.    end if 20.   end for 21.  end for 22.  if $i\ge k⌈e·\sqrt{2\pi \ell }⌉$ then 23.   ${\widehat{K}}^{\left(0\right)}\leftarrow {\widehat{K}}^{\left(best\right)}$ 24.   $k\leftarrow k+1$ 25.  end if 26. end while 27. if$I-{\widehat{K}}_3^{\left(best\right)}{C}_2{B}_2$ is invertible then 28.  return ${\widehat{K}}^{\left(best\right)},f^{\left(best\right)}$ 29. else 30.  return $\mathrm{failure}$ 31. end if

Remark 3.

Please note that ${\widehat{K}}^{\left(0\right)}$ should satisfy the feasibility condition that $I-{\widehat{K}}_3^{\left(0\right)}{C}_2{B}_2$ is invertible, in order to raise the probability that $I-{\widehat{K}}_3^{\left(best\right)}{C}_2{B}_2$ would be invertible. Indeed, one can execute the algorithm without this condition, i.e., execute the algorithm from an infeasible point.

Note also that during the optimization in the main loop of Algorithm 1, we do not check this feasibility condition and we may pass to a better point which might not be feasible, because such points might lead to a better feasible points. This is made because the set of non-feasible points has measure 0 and thus, we should not be disturbed by the possibility of temporarily stepping on infeasible points, during the search. The feasibility condition is checked at the end of Algorithm 1 and experience with the algorithm shows that almost always we end up with a feasible point, when we start from a feasible point (which is almost always the case, if ${\widehat{K}}^{\left(0\right)}$ was found).

5. Experimental Section

To synthesize LQR or $H_{\infty }$ optimal PID controller, we execute the Ray-Shooting algorithm ([6]) in order to get a starting SOF ${\widehat{K}}^{\left(0\right)}$ such that the augmented system (1) is $\Omega$-stable and $I-{\widehat{K}}_3^{\left(0\right)}{C}_2{B}_2$ is invertible. Next, we apply Algorithm 1 in order to get ${\widehat{K}}^{\left(best\right)}$ from which the PID controller is derived, using (2). For the $H_2$ optimal PID controller, we execute the Ray-Shooting algorithm, in order to get a starting SOF ${\widehat{K}}^{\left(0\right)}$, as above. Next, we execute the algorithm from [7], in order to get a SOF ${\widehat{K}}^{\left(1\right)}$ with the structure (6). This step is the most demanding step as can be seen from the following example. Finally, we execute Algorithm 1 in order to get ${\widehat{K}}^{\left(best\right)}$ from which the PID controller is derived, using (2).

Example 1.

For the AC1 system (taken from COMPl${}_e$ib—see [27,28,29]) with the following state-space model:

$$\begin{array}{c}A=\left[\begin{array}{ccccc}\hfill 0& \hfill 0& \hfill 1.132& \hfill 0& \hfill -1\\ \hfill 0& \hfill -0.0538& \hfill -0.1712& \hfill 0& \hfill 0.0705\\ \hfill 0& \hfill 0& \hfill 0& \hfill 1& \hfill 0\\ \hfill 0& \hfill 0.0485& \hfill 0& \hfill -0.8556& \hfill -1.013\\ \hfill 0& \hfill -0.2909& \hfill 0& \hfill 1.0532& \hfill -0.6859\end{array}\right],\hfill \\ B_1=\left[\begin{array}{ccc}\hfill 0.03593& \hfill 0& \hfill 0.01672\\ \hfill 0& \hfill 0.00989& \hfill 0\\ \hfill 0& \hfill -0.07548& \hfill 0\\ \hfill 0& \hfill 0& \hfill 0.05635\\ \hfill 0.00145& \hfill 0& \hfill 0.06743\end{array}\right],B_2=\left[\begin{array}{ccc}\hfill 0& \hfill 0& \hfill 0\\ \hfill -0.12& \hfill 1& \hfill 0\\ \hfill 0& \hfill 0& \hfill 0\\ \hfill 4.419& \hfill 0& \hfill -1.665\\ \hfill 1.575& \hfill 0& \hfill -0.0732\end{array}\right]\hfill \\ C_1=\left[\begin{array}{ccccc}\hfill 0& \hfill 0.707106781186547& \hfill 0& \hfill 0& \hfill 0\\ \hfill 0& \hfill 0& \hfill 0.707106781186547& \hfill 0& \hfill 0\end{array}\right.,\hfill \end{array}$$

and

$$\begin{array}{c}{D}_{1,1}=0_{2\times 3}\hfill \\ D_{1,2}=\left[\begin{array}{ccc}\hfill 0.707106781186547& \hfill 0& \hfill 0\\ \hfill 0& \hfill 0.707106781186547& \hfill 0\end{array}\right]\hfill \\ C_2=\left[\begin{array}{ccccc}\hfill 1& \hfill 0& \hfill 0& \hfill 0& \hfill 0\\ \hfill 0& \hfill 1& \hfill 0& \hfill 0& \hfill 0\\ \hfill 0& \hfill 0& \hfill 1& \hfill 0& \hfill 0\end{array}\right]\hfill \\ D_{2,1}=0_{3\times 3}\hfill \end{array}$$

with

$$\sigma \left(A\right)=\left\{\begin{array}{c}\hfill 0\\ \hfill -0.780052457480417\pm 1.029636749384622i\\ \hfill -0.017597542519583\pm 0.182585375210158i\end{array}\right\},$$

let the objective be to find a PID LQR-optimal controller, where the closed-loop eigenvalues should be in the rectangular region defined by:

$$\Omega =\left\{z\in \mathbb{C}\left|-1\le \Re \left(z\right)\le -0.1,\left|\Im \left(z\right)\right|\le 1\right.\right\}.$$

With the following parameters for the algorithm: $m=100,n=100,s=100,h=r_{\infty }=100,ϵ={10}^{-16},Q=I_8,R=I_2$ we had the following feasible SOF

$${\widehat{K}}^{\left(0\right)}=\left[\begin{array}{ccc}{\widehat{K}}_1^{\left(0\right)}& {\widehat{K}}_2^{\left(0\right)}& {\widehat{K}}_3^{\left(0\right)}\end{array}\right],$$

where

$$\begin{array}{c}{\widehat{K}}_1^{\left(0\right)}=\left[\begin{array}{ccc}\hfill -1.139087350327457& \hfill -0.076805805847702& \hfill -1.175733971595813\\ \hfill 0.175620269295905& \hfill 1.249662743416654& \hfill -0.394592760970558\\ \hfill -2.524143751375688& \hfill -0.262506326222373& \hfill -3.380417685135832\end{array}\right]\hfill \\ {\widehat{K}}_2^{\left(0\right)}=\left[\begin{array}{ccc}\hfill -0.265188982290714& \hfill 0.082432589634872& \hfill -0.236137372902669\\ \hfill 0.106177753653144& \hfill 0.443637620004996& \hfill -0.074457362520092\\ \hfill -0.527328403725350& \hfill 0.200910957443737& \hfill -0.850599462623321\end{array}\right]\hfill \\ {\widehat{K}}_3^{\left(0\right)}=\left[\begin{array}{ccc}\hfill -1.120493844818183& \hfill 0.190581999014316& \hfill -0.522738193473694\\ \hfill -0.036425300944087& \hfill -0.028600461083115& \hfill -0.183279009240151\\ \hfill -3.229623103192055& \hfill 0.550761163041108& \hfill -2.205172558841497\end{array}\right],\hfill \end{array}$$

found by the Ray-Shooting algorithm in $CPU-time=0.0001\left[sec\right]$ and with LQR functional value $f_{LQR}\left({\widehat{K}}^{\left(0\right)}\right)=16.746246825328360$. Executing Algorithm 1, we had

$${\widehat{K}}^{\left(best\right)}=\left[\begin{array}{ccc}{\widehat{K}}_1^{\left(best\right)}& {\widehat{K}}_2^{\left(best\right)}& {\widehat{K}}_3^{\left(best\right)}\end{array}\right],$$

where

$$\begin{array}{c}{\widehat{K}}_1^{\left(best\right)}=\left[\begin{array}{ccc}\hfill -1.085160433924800& \hfill -0.152380406023858& \hfill -1.111840435453025\\ \hfill 0.173224251962231& \hfill 1.253268146513151& \hfill -0.449756208653024\\ \hfill -2.518304745419560& \hfill -0.284285880826183& \hfill -3.437786815141788\end{array}\right]\hfill \\ {\widehat{K}}_2^{\left(best\right)}=\left[\begin{array}{ccc}\hfill -0.221989603556075& \hfill 0.073112868340369& \hfill -0.241805447889791\\ \hfill 0.119158793726812& \hfill 0.403709017174812& \hfill -0.120711030008226\\ \hfill -0.541730211464951& \hfill 0.187380061037601& \hfill -0.903873015737950\end{array}\right]\hfill \\ {\widehat{K}}_3^{\left(best\right)}=\left[\begin{array}{ccc}\hfill -1.068327578247724& \hfill 0.146908957308029& \hfill -0.525373002842116\\ \hfill -0.006608329355743& \hfill 0.015609869717001& \hfill -0.182159600265849\\ \hfill -3.166355434538844& \hfill 0.515059657329335& \hfill -2.156140670934107\end{array}\right],\hfill \end{array}$$

with $f_{LQR}\left({\widehat{K}}^{\left(best\right)}\right)=13.601550793243616$ (a $18.778512\%$ relative improvement), in $CPU-time=76.343750\left[sec\right]$. The resulting PID controller matrices are:

$$\begin{array}{c}{K}_P=\left[\begin{array}{ccc}\hfill -1.040671697032480& \hfill 0.034045803431248& \hfill -1.158219284410088\\ \hfill 0.177951420403750& \hfill 1.273076937581805& \hfill -0.454684211704099\\ \hfill -2.362328186885230& \hfill 0.369320416016725& \hfill -3.600390066086302\end{array}\right]\hfill \\ K_I=\left[\begin{array}{ccc}\hfill -0.200613800205014& \hfill 0.131015509566988& \hfill -0.255249027751546\\ \hfill 0.121430088217771& \hfill 0.409861485359054& \hfill -0.122139482935926\\ \hfill -0.466787102815339& \hfill 0.390385480823640& \hfill -0.951005918753228\end{array}\right]\hfill \\ K_D=\left[\begin{array}{ccc}\hfill -1.050500770135795& \hfill 0.146612915752025& \hfill -0.542836776832826\\ \hfill -0.004714134766140& \hfill 0.015578413703667& \hfill -0.184015220520444\\ \hfill -3.103855027173255& \hfill 0.514021741975605& \hfill -2.217368288447396\end{array}\right]\hfill \end{array}$$

and the resulting closed-loop eigenvalues are

$$\sigma \left({\widehat{A}}_{c\ell }\left({\widehat{K}}^{\left(best\right)}\right)\right)=\left\{\begin{array}{c}\hfill -0.915131493828863\pm 0.914610403731497i\\ \hfill -0.881963764205315\pm 0.191516911259354i\\ \hfill -0.730009436848610\\ \hfill -0.616780364831122\\ \hfill -0.329413066341451\pm 0.071827905231358i\end{array}\right\}.$$

Example 2.

For the same system as in Example 1 with the objective to find a PID $H_{\infty }$-norm optimal controller, where the closed-loop eigenvalues should be in the same rectangular region Ω, starting from the same ${\widehat{K}}^{\left(0\right)}$ with $f_{H_{\infty }}\left({\widehat{K}}^{\left(0\right)}\right)=0.095033139822324$, we had

$${\widehat{K}}^{\left(best\right)}=\left[\begin{array}{ccc}{\widehat{K}}_1^{\left(best\right)}& {\widehat{K}}_2^{\left(best\right)}& {\widehat{K}}_3^{\left(best\right)}\end{array}\right],$$

where

$$\begin{array}{c}{\widehat{K}}_1^{\left(best\right)}=\left[\begin{array}{ccc}\hfill -1.087304300685670& \hfill -0.045931838436646& \hfill -1.085574892719301\\ \hfill 0.133605778924987& \hfill 1.184070254573485& \hfill -0.359343489961314\\ \hfill -2.678924736751655& \hfill -0.264961819769927& \hfill -3.490264530663762\end{array}\right]\hfill \\ {\widehat{K}}_2^{\left(best\right)}=\left[\begin{array}{ccc}\hfill -0.275424882288615& \hfill 0.138303327101369& \hfill -0.075398358659212\\ \hfill 0.090631547438998& \hfill 0.331640807629037& \hfill -0.054696779167659\\ \hfill -0.516085927047396& \hfill 0.184332519865895& \hfill -0.733373655058220\end{array}\right]\hfill \\ {\widehat{K}}_3^{\left(best\right)}=\left[\begin{array}{ccc}\hfill -1.165473798084745& \hfill 0.076978746950025& \hfill -0.595094516639020\\ \hfill -0.052626935909114& \hfill -0.061681173045286& \hfill -0.250098357975546\\ \hfill -3.316496450850004& \hfill 0.694282662038967& \hfill -2.314168285879468\end{array}\right],\hfill \end{array}$$

with $f_{H_{\infty }}\left({\widehat{K}}^{\left(best\right)}\right)=0.072175978563672$ (a $24.051779\%$ relative improvement), in $CPU-time=126\left[sec\right]$. The resulting PID controller matrices are:

$$\begin{array}{c}{K}_P=\left[\begin{array}{ccc}\hfill -1.068321790034687& \hfill 0.039576561177247& \hfill -1.102041006775961\\ \hfill 0.118395560901336& \hfill 1.115554481033627& \hfill -0.346149598520451\\ \hfill -2.507718675999537& \hfill 0.506251036541779& \hfill -3.638774839724999\end{array}\right]\hfill \\ K_I=\left[\begin{array}{ccc}\hfill -0.266534452864325& \hfill 0.160949047585517& \hfill -0.078679654648945\\ \hfill 0.083507865416445& \hfill 0.313495349833722& \hfill -0.052067557547737\\ \hfill -0.435901834009705& \hfill 0.388577609900052& \hfill -0.762968147108973\end{array}\right]\hfill \\ K_D=\left[\begin{array}{ccc}\hfill -1.159203620514524& \hfill 0.071881042424526& \hfill -0.607938725615309\\ \hfill -0.057651074688413& \hfill -0.057596508025015& \hfill -0.239806609825754\\ \hfill -3.259944795394558& \hfill 0.648305713745055& \hfill -2.430012099310619\end{array}\right]\hfill \end{array}$$

and the resulting closed-loop eigenvalues are

$$\sigma \left({\widehat{A}}_{c\ell }\left({\widehat{K}}^{\left(best\right)}\right)\right)=\left\{\begin{array}{c}\hfill -0.965762603991179\pm 0.955842860416991i\\ \hfill -0.507539612520429\pm 0.539780823728534i\\ \hfill -0.929257645094484\pm 0.106820410747074i\\ \hfill -0.402381802562664\pm 0.042181153576293i\end{array}\right\}.$$

Example 3.

For the system from Example 1, with the same Ω and the same ${\widehat{K}}^{\left(0\right)}$, the algorithm from [7] failed to find structured SOF ${\widehat{K}}^{\left(1\right)}$. Indeed, from another starting SOF, the above-mentioned algorithm succeeded to find a structured SOF. For

$${\widehat{K}}^{\left(0\right)}=\left[\begin{array}{ccc}{\widehat{K}}_1^{\left(0\right)}& {\widehat{K}}_2^{\left(0\right)}& {\widehat{K}}_3^{\left(0\right)}\end{array}\right],$$

where

$$\begin{array}{c}{\widehat{K}}_1^{\left(0\right)}=\left[\begin{array}{ccc}\hfill -0.001662379279286& \hfill -0.330320315815391& \hfill -0.204742383931633\\ \hfill 0.107146675564809& \hfill 0.994156785847185& \hfill 0.016601261567192\\ \hfill -0.352905474519842& \hfill -0.942806073864177& \hfill 0.096372856789800\end{array}\right]\hfill \\ {\widehat{K}}_2^{\left(0\right)}=\left[\begin{array}{ccc}\hfill 0.021994917878383& \hfill -0.051106382240571& \hfill 0.041069916425995\\ \hfill 0.041186023579647& \hfill 0.283840911655084& \hfill 0.016793398034650\\ \hfill 0.000278730967043& \hfill -0.136566568491544& \hfill 0.140207270149126\end{array}\right]\hfill \\ {\widehat{K}}_3^{\left(0\right)}=\left[\begin{array}{ccc}\hfill -0.066971504986511& \hfill 0.041204949529465& \hfill 0.574933673498160\\ \hfill -0.152954829804652& \hfill -0.044219575277023& \hfill 0.285308855216965\\ \hfill -1.286404253117917& \hfill 0.132606565745259& \hfill 1.514386124380330\end{array}\right],\hfill \end{array}$$

found by the Ray-Shooting algorithm from [6], in $CPU-time=0.0001\left[sec\right]$, applying the algorithm from [7] we have found the structured SOF

$${\widehat{K}}^{\left(1\right)}=\left[\begin{array}{ccc}{\widehat{K}}_1^{\left(1\right)}& {\widehat{K}}_2^{\left(1\right)}& {\widehat{K}}_3^{\left(1\right)}\end{array}\right],$$

where

$$\begin{array}{c}{\widehat{K}}_1^{\left(1\right)}=\left[\begin{array}{ccc}\hfill -0.001662379279286& \hfill -0.330320315815391& \hfill -0.204742383931633\\ \hfill 0.107146675564809& \hfill 0.994156785847185& \hfill 0.016601261567192\\ \hfill -0.352905474519842& \hfill -0.942806073864177& \hfill 0.096372856789800\end{array}\right]\hfill \\ {\widehat{K}}_2^{\left(1\right)}=\left[\begin{array}{ccc}\hfill 0.021994917878383& \hfill -0.051106382240571& \hfill 0.041069916425995\\ \hfill 0.041186023579647& \hfill 0.283840911655084& \hfill 0.016793398034650\\ \hfill 0.000278730967043& \hfill -0.136566568491544& \hfill 0.140207270149126\end{array}\right]\hfill \\ {\widehat{K}}_3^{\left(1\right)}=\left[\begin{array}{ccc}\hfill -0.000000000000001& \hfill 1.630868948596140& \hfill 0.213689638336193\\ \hfill -0.000000000000003& \hfill 0.516222895350342& \hfill 0.067639698397123\\ \hfill -3.828454587894386& \hfill 1.446059480160095& \hfill 0.287354792999600\end{array}\right],\hfill \end{array}$$

in $CPU-time=2126.84375\left[sec\right]=35.447395\left[min\right]$, with objective-function value $f_{H_2}\left({\widehat{K}}^{\left(1\right)}\right)=0.364564164927770$ (note the small values of the $1,1$ and the $2,1$ entries of ${\widehat{K}}_3^{\left(1\right)}$ obtained from the structural constraints). Finally, executing Algorithm 1 with ${\widehat{K}}^{\left(1\right)}$ as a starting SOF, we have found

$${\widehat{K}}^{\left(best\right)}=\left[\begin{array}{ccc}{\widehat{K}}_1^{\left(best\right)}& {\widehat{K}}_2^{\left(best\right)}& {\widehat{K}}_3^{\left(best\right)}\end{array}\right],$$

where

$$\begin{array}{c}{\widehat{K}}_1^{\left(best\right)}=\left[\begin{array}{ccc}\hfill 0.076825955660614& \hfill -0.288411611380877& \hfill -0.080218436303677\\ \hfill 0.092501965527188& \hfill 1.046837016022694& \hfill 0.126682619794964\\ \hfill -0.341825773843233& \hfill -0.939476093514249& \hfill -0.014192822579841\end{array}\right]\hfill \\ {\widehat{K}}_2^{\left(best\right)}=\left[\begin{array}{ccc}\hfill 0.054840679045778& \hfill -0.080221174899522& \hfill 0.114230924697635\\ \hfill 0.052402918799798& \hfill 0.301314353311879& \hfill -0.093055199101202\\ \hfill -0.061395620772426& \hfill -0.248330164720565& \hfill 0.198088337294925\end{array}\right]\hfill \\ {\widehat{K}}_3^{\left(best\right)}=\left[\begin{array}{ccc}\hfill -0.000000000000001& \hfill 1.623084863796692& \hfill 0.212669704596573\\ \hfill -0.000000000000003& \hfill 0.692852565144541& \hfill 0.090783146121880\\ \hfill -3.950508497484917& \hfill 1.418275521031274& \hfill -0.019182452364697\end{array}\right],\hfill \end{array}$$

with $f_{H_2}\left({\widehat{K}}^{\left(best\right)}\right)=0.212807638848134$ (a $41.626835\%$ relative improvement), in $CPU-time=99.890625\left[sec\right]$. The resulting PID controller matrices are:

$$\begin{array}{c}{K}_P=\left[\begin{array}{ccc}\hfill 0.346143328950049& \hfill 3.208734600167995& \hfill 0.360572630013263\\ \hfill 0.207466522885631& \hfill 2.539677456718627& \hfill 0.314844821170003\\ \hfill -0.106492276227244& \hfill 2.116381906236221& \hfill 0.370976920044842\end{array}\right]\hfill \\ K_I=\left[\begin{array}{ccc}\hfill 0.203018492396510& \hfill 0.925289371971877& \hfill -0.231015500114683\\ \hfill 0.115656159814819& \hfill 0.730540561515251& \hfill -0.240431886773992\\ \hfill 0.068084340579495& \hfill 0.630299798860383& \hfill -0.103593086737016\end{array}\right]\hfill \\ K_D=\left[\begin{array}{ccc}\hfill -0.000000000000010& \hfill 3.233767462898737& \hfill 0.423714364176847\\ \hfill -0.000000000000007& \hfill 1.380410927195360& \hfill 0.180872602940673\\ \hfill -3.950508497484925& \hfill 2.825713760036134& \hfill 0.165231487536706\end{array}\right]\hfill \end{array}$$

and the resulting closed-loop eigenvalues are

$$\sigma \left({\widehat{A}}_{c\ell }\left({\widehat{K}}^{\left(best\right)}\right)\right)=\left\{\begin{array}{c}\hfill -0.350618286842131\pm 0.990354187385021i\\ \hfill -0.963555748849909\pm 0.645732021985088i\\ \hfill -0.127768318435090\pm 0.216451742947409i\\ \hfill -0.159486067934620\\ \hfill -0.462031919768989\end{array}\right\}.$$

Please note that

$$\begin{array}{cc}\hfill {\widehat{D}}_{1,1}& =D_{1,1}-D_{1,2}{\widehat{K}}_3^{\left(1\right)}{C}_2{B}_1\hfill \\ & ={10}^{-16}·\left[\begin{array}{ccc}\hfill 0.254063466480327& \hfill 0.173472347597681& \hfill 0.118228253814391\\ \hfill 0.762190399440980& \hfill -0.273218947466347& \hfill 0.354684761443172\end{array}\right]\hfill \end{array}$$

and

$$\begin{array}{cc}\hfill {\widehat{D}}_{1,1}& =D_{1,1}-D_{1,2}{\widehat{K}}_3^{\left(best\right)}{C}_2{B}_1\hfill \\ & ={10}^{-16}·\left[\begin{array}{ccc}\hfill 0.254063466480326& \hfill 0.156125112837913& \hfill 0.118228253814391\\ \hfill 0.762190399440979& \hfill -0.277555756156289& \hfill 0.354684761443172\end{array}\right]\hfill \end{array}$$

which we treat numerically as 0.

6. Concluding Remarks

In this article, a randomized method for synthesizing PID optimal controllers (for the standard LQR, $H_{\infty }$ and $H_2$ optimization problems) for coupled MIMO systems was introduced. The method is based on a (direct) reduction of the problem to the problem of SOF (occasionally with some structure). The randomization circumvents the need to solve any SDP or BMI programs. Therefore, the method can be applied to more involved PID controller design that can be formulated as structured SOF problems, e.g., the design of 2-DOF-PID controllers and robust and gain-scheduled PID controllers for Linear Parameter Varying (LPV) systems (see [13]), without the need to solve any BMI programs and thus can be applied to larger systems. Moreover, the method has a proof of convergence (in probability) to a global optimum of the objective function.