Method of Reduction of Variables for Bilinear Matrix Inequality Problems in System and Control Designs

TL;DR

This paper introduces a variable reduction method for bilinear matrix inequality problems in system and control design, transforming complex problems into simpler unconstrained forms to improve solution efficiency.

Contribution

A novel reduction of variables method combining classification, problem transformation, and hybrid search algorithms for BMI problems in control systems.

Findings

Outperforms existing BMI solution methods on benchmark problems.

Effectively addresses feasibility and multiobjective problems.

Reduces decision variables, simplifying complex BMI problems.

Abstract

Bilinear matrix inequality (BMI) problems in system and control designs are investigated in this paper. A solution method of reduction of variables (MRVs) is proposed. This method consists of a principle of variable classification, a procedure for problem transformation, and a hybrid algorithm that combines deterministic and stochastic search engines. The classification principle is used to classify the decision variables of a BMI problem into two categories: 1) external and 2) internal variables. Theoretical analysis is performed to show that when the classification principle is applicable, a BMI problem can be transformed into an unconstrained optimization problem that has fewer decision variables. Stochastic search and deterministic search are then applied to determine the decision variables of the unconstrained problem externally and explore the internal problem structure, respectively. The proposed method can address feasibility, single-objective, and multiobjective problems constrained by BMIs in a unified manner. A number of numerical examples in system and control designs are provided to validate the proposed methodology. Simulations show that the MRVs can outperform existing BMI solution methods in most benchmark problems and achieve similar levels of performance in the remaining problems.

Tables

Table 1. TABLE I: Feasibility Problems

Problems	Existing Solution Methods	Results of MRV
Name		SR $\%$
ST [11]	AM	100
SIP [11]	AM	100
SAFS-I [7]	ILMI	100
SAFS-II [8]	Diffeomorphic state transformations	100
OCS [9]	Two-step procedure	100

Table 2. TABLE II: Single-objective Optimization Problems

Problems	Methods	Results of Existing Methods	Results of MRV
Name			Optimum	Mean	Std	SR $\%$
LPVS (maximization) [29]	BB methods	4.75	4.7575	4.7209	0.0280	100
SSS (maximization) [2]	Path-following methods	1.05	4.1765	3.3347	0.3758	100
MCD (minimization) [15]	Path-following methods	0.7489	0.7600	0.8291	0.0407	100

Table 3. TABLE III: Spectral Abscissa Optimization

Problems		Results of Existing Solution Methods, ${\alpha }_o\left({𝑨}_{𝑭}\right)$							Results of MRV, ${\alpha }_o\left({𝑨}_{𝑭}\right)$
Name	${\alpha }_o\left(𝑨\right)$	HIFOO	LMIRank	PENBMI	CCDM	ICAM	Min	Mean	Min	Mean	Std	SR $\%$
AC1	0.000	-0.2061	-8.4766	-7.0758	-0.8535	-0.7814	-8.4766	-3.4786	-18.0761	-11.8993	3.2210	100
AC4	2.579	-0.0500	-0.0500	-0.0500	-0.0500	-0.0500	-0.05	-0.05	-0.05	-0.05	6.9e-17	100
AC5	0.999	-0.7746	-1.8001	-2.0438	-0.7389	-0.7389	-2.0438	-1.2192	-2.4051	-2.1444	0.1754	100
AC7	0.172	-0.0322	-0.0204	0.0896	-0.0673	-0.0502	-0.0673	-0.0161	-0.0747	-0.0494	0.0088	100
AC8	0.012	-0.1968	-0.4447	0.4447	-0.0755	-0.0640	-0.4447	-0.0672	-0.4447	-0.4447	2.7e-16	100
AC9	0.012	-0.3389	-0.5230	-0.4450	-0.3256	-0.3926	-0.523	-0.405	-2.0823	-0.5776	0.2970	100
AC11	5.451	-0.0003	-5.0577	x	-3.0244	-3.1573	-5.0577	-2.8099	-16.9018	-10.6947	2.6689	100
AC12	0.580	-10.8645	-9.9658	-1.8757	-0.3414	-0.2948	-10.8645	-4.6684	-18.3236	-13.3959	2.8633	100
HE1	0.276	-0.2457	-0.2071	-0.2468	-0.2202	-0.2134	-0.2468	-0.2266	-0.2446	-0.2338	0.0107	100
HE3	0.087	-0.4621	-2.3009	-0.4063	-0.8702	-0.8380	-2.3009	-0.9755	-1.7847	-0.8908	0.3055	100
HE4	0.234	-0.7446	-1.9221	-0.0909	-0.8647	-0.8375	-1.9221	-0.8919	-3.0567	-1.2306	0.5201	100
HE5	0.234	-0.1823	x	-0.2932	-0.0587	-0.0609	-0.2932	-0.1487	-1.1953	-0.6939	0.2193	100
HE6	0.234	-0.0050	-0.0050	-0.0050	-0.0050	-0.0050	-0.005	-0.005	-0.005	-0.005	2.6e-18	100
REA1	1.991	-16.3918	-5.9736	-1.7984	-3.8599	-2.8932	-16.3918	-6.1833	-19.3041	-15.4064	2.4190	100
REA2	2.011	-7.0152	-10.0292	-3.5928	-2.1778	-1.9514	-10.0292	-4.9532	-19.4238	-13.0948	4.2323	100
REA3	0.000	-0.0207	-0.0207	-0.0207	-0.0207	-0.0207	-0.0207	-0.0207	-0.0207	-0.0207	3.5e-15	100
DIS2	1.675	-6.8510	-10.1207	-8.3289	-8.4540	-8.3419	-10.1207	-8.4193	-19.4340	-16.6852	2.5153	100
DIS4	1.442	-36.7203	-0.5420	-92.2842	-8.0989	-5.4467	-92.2842	-28.6184	-16.0222	-11.4094	2.4090	100
WEC1	0.008	-8.9927	-8.7350	-0.9657	-0.8779	-0.8568	-8.9927	-4.0856	-11.9629	-6.1804	2.2291	100
IH	0.000	-0.5000	-0.5000	-0.5000	-0.5000	-0.5000	-0.5	-0.5	-0.1576	-0.0617	0.0407	76.47
CSE1	0.000	-0.4509	-0.4844	-0.4490	-0.2360	-0.2949	-0.4844	-0.383	-0.3489	-0.2282	0.0452	100
TF1	0.000	x	x	-0.0618	-0.1544	-0.0704	-0.1544	-0.0955	-0.2688	-0.1769	0.0396	100
TF2	0.000	x	x	-1.0e-5	-1.0e-5	-1.0e-5	-1.0e-5	-1.0e-5	-1.0e-5	-1.0e-5	1.7e-20	100
TF3	0.000	x	x	-0.0032	-0.0031	-0.0032	-0.0032	-0.0031	-0.0032	-0.0032	8.0e-6	100
NN1	3.606	-3.0458	-4.4021	-4.3358	-0.8746	0.1769	-4.4021	-2.4962	-5.89	-5.6847	0.1812	100
NN5	0.420	-0.0942	-0.0057	-0.0942	-0.0913	-0.0490	-0.0942	-0.0668	-0.094	-0.0915	0.0018	100
NN9	3.281	-2.0789	-0.7048	x	-0.0279	0.0991	-2.0789	-0.6781	-17.8516	-12.1047	2.7270	100
NN13	1.945	-3.2513	-4.5310	-9.0741	-3.4318	-0.2783	-9.0741	-4.1133	-13.6061	-8.5606	4.7341	100
NN15	0.000	-6.9983	-11.0743	-0.0278	-0.8353	-1.0409	-11.0743	-3.9953	-10.9821	-10.3002	0.8034	100
NN17	1.170	-0.6110	-0.5130	x	-0.6008	-0.5991	-0.611	-0.3244	-0.6107	-0.6007	0.0196	100

Table 4. TABLE IV: H 2 subscript 𝐻 2 H_{2} Optimization

Problems		Results of Existing Solution Methods, ${\Vert {𝑮}_{c\ell }\left(𝑭\right)\Vert }_2$					Results of MRV, ${\Vert {𝑮}_{c\ell }\left(𝑭\right)\Vert }_2$
Name	${\Vert {𝑮}_{o\ell }\Vert }_2$	HIFOO	PENBMI	CCDM	Min	Mean	Min	Mean	Std	SR $\%$
AC1	Inf	0.025	0.0061	0.054	0.0061	0.0283	0.015	0.0187	0.0019	100
AC2	Inf	0.0257	0.0075	0.054	0.0075	0.029	0.01566	0.0188	0.0017	100
AC3	25.5798	2.0964	2.0823	2.1117	2.0823	2.0968	2.1206	2.231	0.0836	100
AC4	Inf	11.0269	x	11.0269	11.0269	11.0269	11.0269	11.0269	2.82e-15	100
AC6	24.6067	2.8648	2.8648	2.8664	2.8648	2.8653	3.026	3.6263	0.6533	100
AC7	Inf	0.0172	0.0162	0.0176	0.0162	0.017	0.0162	0.0164	0.0001	100
AC8	Inf	0.633	0.7403	0.6395	0.633	0.6709	0.6321	0.6813	0.0417	100
AC12	Inf	0.0022	0.0106	0.0992	0.0022	0.0373	0.0627	0.1129	0.0295	100
AC15	176.4515	1.5458	1.4811	1.5181	1.4811	1.515	1.6564	1.828	0.1709	100
AC16	176.4515	1.4769	1.4016	1.4427	1.4016	1.4404	1.4641	1.5307	0.0335	100
AC17	10.2650	1.5364	1.5347	1.5507	1.5347	1.5406	1.5392	1.5429	0.0041	100
HE2	13.8541	3.4362	3.4362	4.7406	3.4362	3.871	3.7494	6.1145	1.1487	100
HE3	Inf	0.0197	0.0071	0.1596	0.0071	0.0621	0.0333	0.1026	0.0808	100
HE4	Inf	6.6436	6.5785	7.1242	6.5785	6.7821	15.7738	27.0193	8.7713	100
REA1	Inf	0.9442	0.9422	1.0622	0.9422	0.9828	0.9593	0.9864	0.0176	100
REA2	Inf	1.0339	1.0229	1.1989	1.0229	1.0852	1.0261	1.0319	0.0134	100
DIS1	5.1491	0.6705	0.1174	0.7427	0.1174	0.5102	0.51	0.7455	0.1531	100
DIS2	Inf	0.4013	0.37	0.3819	0.37	0.3844	0.372	0.381	0.0128	100
DIS3	11.6538	0.9527	0.9434	1.0322	0.9434	0.9761	0.997	1.0623	0.0288	100
DIS4	Inf	1.0117	0.9696	1.0276	0.9696	1.0029	1.0644	1.1091	0.0351	100
WEC1	Inf	7.394	8.1032	12.9093	7.394	9.4688	12.1017	16.2366	1.9181	100
WEC2	66.5622	6.7908	7.6502	12.2102	6.7908	8.8837	13.2889	16.5298	1.3581	100
AGS	7.0412	6.9737	6.9737	6.9838	6.9737	6.977	7.1807	10.1753	2.4965	100
BDT1	0.0397	0.0024	x	0.0017	0.0017	0.002	3.52e-05	5.44e-05	1.22e-05	100
MFP	12.6469	6.9724	6.9724	7.0354	6.9724	6.9934	7.0556	7.6688	0.7976	100
PSM	3.8474	0.033	0.0007	0.1753	0.0007	0.0697	0.0217	0.04	0.0149	100
EB2	4.0000	0.064	0.0084	0.1604	0.0084	0.0776	0.0832	0.086	0.0091	100
EB3	1.26e03	0.0732	0.0072	0.0079	0.0072	0.0294	0.0846	0.0918	0.0141	100
TF1	Inf	0.0945	x	0.1599	0.0945	0.1272	0.1949	0.6965	1.2484	100
TF2	Inf	11.1803	x	11.1803	11.1803	11.1803	11.1803	11.1803	1.48e-14	100
TF3	Inf	0.1943	0.1424	0.2565	0.1424	0.1977	0.2568	2.0745	1.6401	97.67
NN2	Inf	1.1892	1.1892	1.1892	1.1892	1.1892	1.1892	1.1892	3.82e-06	100
NN4	5.5634	1.8341	1.8335	1.859	1.8335	1.8422	1.8945	1.989	0.056	100
NN8	5.9220	1.5152	1.5117	1.5725	1.5117	1.5331	1.5241	1.5518	0.017	100
NN11	0.1420	0.1178	0.079	0.1263	0.079	0.1077	0.0972	0.1137	0.0102	100
NN13	Inf	26.1012	26.1314	62.3995	26.1012	38.2107	30.1629	34.4666	3.9562	100
NN14	Inf	26.1448	26.1314	62.3995	26.1314	38.2252	29.6438	35.6657	7.4852	100
NN15	Inf	0.0245	x	0.021	0.021	0.0227	0.0034	0.0035	8.82e-05	100
NN16	Inf	0.1195	0.1195	0.1195	0.1195	0.1195	0.1208	0.2085	0.068	100
NN17	Inf	3.253	3.2404	3.3329	3.2404	3.2754	3.2554	3.2881	0.1843	100

Table 5. TABLE V: H ∞ subscript 𝐻 H_{\infty} Optimization

Problems		Results of Existing Solution Methods, ${\Vert {𝑮}_{c\ell }\left(𝑭\right)\Vert }_{\infty }$					Results of MRV, ${\Vert {𝑮}_{c\ell }\left(𝑭\right)\Vert }_{\infty }$
Name	${\Vert {𝑮}_{o\ell }\Vert }_{\infty }$	HIFOO	PENBMI	CCDM	Min	Mean	Min	Mean	Std	SR $\%$
AC1	2.1672	0.0000	x	0.0177	0.0000	0.0088	0.0405	0.0907	0.0285	100
AC2	2.1672	0.1115	x	0.1140	0.1115	0.1127	0.1262	0.1917	0.0310	100
AC3	352.6869	4.7021	x	3.4859	3.4859	4.094	3.9206	4.5709	0.4217	100
AC4	69.9900	0.9355	x	69.9900	0.9355	35.4627	69.99	69.99	1.29e-13	100
AC6	391.7820	4.1140	x	4.1954	4.114	4.1547	4.8138	6.9232	3.1554	100
AC7	0.0424	0.0651	0.3810	0.0548	0.0548	0.1669	0.0315	0.0316	6.27e-05	100
AC8	1.7e03	2.0050	x	3.0520	2.005	2.5285	1.4305	1.8223	0.4017	100
AC9	Inf	1.0048	x	0.9237	0.9237	0.9642	3.2926	5.1355	1.069	100
AC11	Inf	3.5603	x	3.0104	3.0104	3.28535	3.1158	4.0119	0.5472	100
AC12	586.9176	0.3160	x	2.3025	0.316	1.3092	1.3532	1.9379	0.1729	100
AC15	2.4e03	15.2074	427.4106	15.1995	15.1995	152.6058	17.1925	18.2818	0.4480	100
AC16	2.4e03	15.4969	x	14.9881	14.9881	15.2425	15.8600	16.6389	0.5547	100
AC17	30.8328	6.6124	x	6.6373	6.6124	6.6248	6.6124	6.6124	1.30e-06	100
HE1	0.5598	0.1540	1.5258	0.1807	0.154	0.6201	0.1538	0.1595	0.0045	100
HE2	81.8318	4.4931	x	6.7846	4.4931	5.6388	4.3681	5.5034	0.7626	100
HE3	1.4618	0.8545	1.6843	0.9243	0.8545	1.1543	0.8570	0.9142	0.0381	100
HE4	174.2975	23.3448	x	22.8713	22.8713	23.108	46.5677	65.3844	7.8214	100
HE5	2.0802	8.8952	x	37.3906	8.8952	23.1429	20.8784	137.7817	155.8509	100
REA1	25.7708	0.8975	x	0.8815	0.8815	0.8895	0.8836	0.9073	0.0309	100
REA2	26.3449	1.1881	x	1.4188	1.1881	1.3034	1.1471	1.168	0.0125	100
REA3	Inf	74.2513	74.446	74.5478	74.2513	74.415	74.2513	75.5692	2.3953	100
DIS1	17.3209	4.1716	x	4.1943	4.1716	4.1829	4.3197	4.7678	0.4376	100
DIS2	0.9016	1.0548	1.7423	1.1546	1.0548	1.3172	1.0604	1.1364	0.039	100
DIS3	32.0698	1.0816	x	1.1382	1.0816	1.1099	1.2727	1.3733	0.0436	100
DIS4	3.1304	0.7465	x	0.7498	0.7465	0.7481	0.9486	1.0203	0.0411	100
TG1	130.3418	12.8462	x	12.9336	12.8462	12.8899	14.2157	25.1589	13.8431	100
AGS	8.1820	8.1732	188.0315	8.1732	8.1732	68.126	10.0239	20.963	6.2165	100
WEC2	354.3162	4.2726	32.9935	6.6082	4.2726	14.6247	7.8382	10.3568	1.5282	100
WEC3	180.0408	4.4497	200.1467	6.8402	4.4497	70.4788	7.2021	9.6854	1.4185	100
BDT1	5.1426	0.2664	x	0.8562	0.2664	0.5613	0.2662	0.2669	0.0008	100
MFP	83.1407	31.5899	x	31.6079	31.5899	31.5989	33.9193	51.8236	25.4745	100
IH	Inf	1.9797	x	1.1858	1.1858	1.5827	30.1004	450.2228	1231.8531	90.38
CSE1	1.3e13	0.0201	x	0.0220	0.0201	0.021	0.0198	0.0199	2.00e-05	100
PSM	4.2328	0.9202	x	0.9227	0.9202	0.9214	0.9202	0.9208	0.001	100
EB1	39.9526	3.1225	39.9526	2.0276	2.0276	15.0342	1.888	1.888	4.81e-08	100
EB2	39.9526	2.0201	39.9547	0.8148	0.8148	14.2632	0.8142	0.8142	8.44e-16	100
EB3	3.9e06	2.0575	3995311.074	0.8153	0.8153	1331771.316	0.8143	0.8143	6.17e-16	100
NN1	Inf	13.9782	14.6882	18.4813	13.9782	15.7159	15.5294	16.6317	0.8991	100
NN2	Inf	2.2216	x	2.2216	2.2216	2.2216	2.2038	2.2056	0.0021	100
NN4	31.0435	1.3627	x	1.3802	1.3627	1.3714	1.4327	1.6037	0.0855	100
NN8	46.5086	2.8871	78281181.15	2.9345	2.8871	26093728.99	2.9193	2.9977	0.047	100
NN9	3.7675	28.9083	x	32.1222	28.9083	30.5152	30.7173	35.299	7.047	100
NN11	0.1703	0.1037	x	0.1566	0.1037	0.1301	0.1075	0.1374	0.0127	100
NN15	Inf	0.1039	x	0.1194	0.1039	0.1116	0.098	0.0982	0.0001	100
NN16	6.4e14	0.9557	x	0.9656	0.9557	0.9606	2.3044	6.8293	3.2719	100
NN17	2.8284	11.2182	x	11.2381	11.2182	11.2281	11.2042	11.6262	0.5366	100

Table 6. TABLE VI: Additional Problems in Spectral Abscissa, H 2 subscript 𝐻 2 H_{2} , and H ∞ subscript 𝐻 H_{\infty} Optimization

Problems		Results of MRV, ${\alpha }_o\left({𝑨}_{𝑭}\right)$
Name	${\alpha }_o\left(𝑨\right)$	Min	Mean	Std	SR $\%$
AC18	0.1015	-1.9248	-1.1526	0.3836	100
DIS5	1.0192	-2.7044	-2.3709	0.2454	100
PAS	0	-2.05e-05	-1.61e-05	5.87e-06	100
NN12	1.0000	-2.4761	-1.9860	0.5362	100
Problems		Results of MRV, ${\Vert {𝑮}_{c\ell }\left(𝑭\right)\Vert }_2$
Name	${\Vert {𝑮}_{o\ell }\Vert }_2$	Min	Mean	Std	SR $\%$
AC18	Inf	20.0248	21.1601	0.6287	100
DIS5	Inf	0.0013	0.0019	0.0006	100
NN12	Inf	8.6989	10.5373	3.7473	100
Problems		Results of MRV, ${\Vert {𝑮}_{c\ell }\left(𝑭\right)\Vert }_{\infty }$
Name	${\Vert {𝑮}_{o\ell }\Vert }_{\infty }$	Min	Mean	Std	SR $\%$
AC18	140.3365	10.8088	11.9210	1.7586	100
DIS5	0.0108	28.7928	29.3512	0.4921	100
NN12	Inf	22.4556	40.5618	28.6066	100

Equations

$$\left[\begin{array}[]{cc}\bm{A}&\bm{B}^{T}\\ \bm{B}&\bm{C}\\ \end{array}\right]=\left[\begin{array}[]{cc}\bm{A}&\star\\ \bm{B}&\bm{C}\\ \end{array}\right]=\left[\begin{array}[]{cc}\bm{A}&\bm{B}^{T}\\ \star&\bm{C}\\ \end{array}\right].$$

$$\begin{array}{rl}\underset{\mathbf{\omega }}{\min }& \mathbf{f}(\mathbf{\omega })\\ \text{mbox}subjectto& \mathbf{\omega }\in \Omega \end{array}$$

$${\mathcal{P}}^{\ast }=\{\mathbf{\omega }\in \Omega :\nexists {\mathbf{\omega }}^{\prime }\in \Omega \text{mbox}suchthat{\mathbf{\omega }}^{\prime }{⪯}_{\mathbf{f}}\mathbf{\omega }\}.$$

$$\mathcal{BM}\mathcal{I}(\mathbf{\alpha },\mathbf{X})<0$$

$$\begin{array}{rl}\underset{\mathbf{\alpha },\mathbf{X}}{\min }& \mathcal{F}(\mathbf{\alpha })\\ \text{mbox}subjectto& \mathcal{BM}\mathcal{I}(\mathbf{\alpha },\mathbf{X})<0\end{array}$$

$$\dot{\mathbf{x}}(t)=\sum _{i=1}^2{\xi }_i(\mathbf{x}(t)){\mathbf{A}}_i\mathbf{x}(t)+\mathbf{p},{\mathbf{p}}^T\mathbf{p}\le {\mu }^2\mathbf{x}(t{)}^T\mathbf{x}(t).$$

$${\mathbf{A}}_{\ell }^T{\mathbf{P}}_i+{\mathbf{P}}_i{\mathbf{A}}_{\ell }+{\mu }^2\mathbf{I}-\sum _{j=1}^2{\tau }_{\ell ij}({\mathbf{P}}_j-{\mathbf{P}}_i)<0,\text{mbox}for\ell ,i=1,2$$

$$\dot{\mathbf{x}}(t)=\mathbf{A}(t)\mathbf{x}(t),\mathbf{A}(t)\in \text{mbox}convexhull\{{\mathbf{A}}_1,{\mathbf{A}}_2\}$$

$$\bm{A}_{1}=\left[\begin{array}[]{cc}0&1\\ -2&-1\\ \end{array}\right]\mbox{ and }\bm{A}_{2}=\left[\begin{array}[]{cc}0&1\\ -2-\varsigma&-1\\ \end{array}\right].$$

$$\begin{array}{rl}& (1-{\delta }_2)({\mathbf{P}}_2{\mathbf{A}}_1,\star )+{\delta }_2({\mathbf{P}}_2-{\mathbf{P}}_1)<0\\ & (1-{\delta }_1)({\mathbf{P}}_1{\mathbf{A}}_2,\star )-{\delta }_1({\mathbf{P}}_2-{\mathbf{P}}_1)<0\\ & ({\mathbf{P}}_1{\mathbf{A}}_1,\star )<0,({\mathbf{P}}_2{\mathbf{A}}_2,\star )<0\\ & 0<{\mathbf{P}}_i<\mathbf{I},0\le {\delta }_i\le 1,i=1,2.\end{array}$$

$$\begin{array}{rl}\underset{\varsigma ,{\delta }_i,{\mathbf{P}}_i}{\max }& \varsigma \\ \text{mbox}subjectto& (\text{ref}e{q}_{L}PV{S}_{B}MI).\end{array}$$

$$\begin{array}{rl}\underset{\beta ,\mathbf{F},\mathbf{P}}{\min }& -\sigma \beta +\sum _i\sum _j|[\mathbf{F}{]}_{ij}|\\ \text{mbox}subjectto& (\mathbf{P}{\mathbf{A}}_{\mathbf{F}},\star )+2\beta \mathbf{P}<0,\mathbf{P}>0\end{array}$$

$$\begin{split}\min_{\beta,\bm{F},\bm{P}}&\;\left[\begin{array}[]{cc}-\beta&\sum_{i}\sum_{j}|[\bm{F}]_{ij}|\\ \end{array}\right]^{T}\\ \mbox{subject to }&\;(\bm{P}\bm{A}_{\bm{F}},\star)+2\beta\bm{P}<0,\bm{P}>0.\end{split}$$

$$(\lambda^{*}(\bm{\alpha}),\bm{X}^{*}(\bm{\alpha}))=\begin{array}[]{l}\arg_{\lambda,\bm{X}}\min_{\lambda,\bm{X}}\;\lambda\\ \mbox{subject to }\mathcal{BMI}(\bm{\alpha},\bm{X})<\lambda\bm{I}.\end{array}$$

$$\begin{array}{rl}\underset{\mathbf{\alpha }}{\min }& \mathcal{F}(\mathbf{\alpha })\\ \text{mbox}subjectto& {\lambda }^{\ast }(\mathbf{\alpha })<0.\end{array}$$

$${\mathbf{\alpha }}^{\prime }{⪯}_{\mathcal{F}}\tilde{\mathbf{\alpha }},{\lambda }^{\ast }(\tilde{\mathbf{\alpha }})<0,\text{mbox}and{\lambda }^{\ast }({\mathbf{\alpha }}^{\prime })<0.$$

$$({\mathbf{\alpha }}^{\prime },{\mathbf{X}}^{\ast }({\mathbf{\alpha }}^{\prime })){⪯}_{\mathcal{F}}(\tilde{\mathbf{\alpha }},\tilde{\mathbf{X}})$$

$$({\mathbf{\alpha }}^{\prime },{\mathbf{X}}^{\prime }){⪯}_{\mathcal{F}}(\tilde{\mathbf{\alpha }},\tilde{\mathbf{X}}),\mathcal{BM}\mathcal{I}({\mathbf{\alpha }}^{\prime },{\mathbf{X}}^{\prime })<0,\text{mbox}and\mathcal{BM}\mathcal{I}(\tilde{\mathbf{\alpha }},\tilde{\mathbf{X}})<0.$$

$$\tilde{\bm{\mathcal{F}}}(\bm{\alpha})=\left[\begin{array}[]{cc}\bm{\mathcal{F}}(\bm{\alpha})^{T}&\max\{0,\lambda^{*}(\bm{\alpha})\}\\ \end{array}\right]^{T}$$

$$\underset{\mathbf{\alpha }}{\min }\tilde{\mathcal{F}}(\mathbf{\alpha }).$$

$$h(\mathbf{q},{\mathbf{\lambda }}^{pre})=\frac{1}{2}||eig(\mathbf{A}(\mathbf{q}))-{\mathbf{\lambda }}^{pre}|{|}_2^2$$

$${\mathbf{q}}^{\ast }({\mathbf{\lambda }}^{pre})={\arg }_{\mathbf{q}}\underset{\mathbf{q}}{\min }h(\mathbf{q},{\mathbf{\lambda }}^{pre})$$

$$\begin{array}{rl}& \mathbf{A}(\mathbf{q}){\mathbf{X}}_e(\mathbf{q})\\ =& {\mathbf{X}}_e(\mathbf{q})diag([eig(\mathbf{A}(\mathbf{q})){]}_1,[eig(\mathbf{A}(\mathbf{q})){]}_2,...,[eig(\mathbf{A}(\mathbf{q})){]}_{n_x}).\end{array}$$

$$\frac{\partial[eig(\bm{A}(\bm{q}))]_{i}}{\partial[\bm{q}]_{m}}=\left[\begin{array}[]{c}\bm{X}_{e}(\bm{q})^{-1}\frac{\partial\bm{A}(\bm{q})}{\partial[\bm{q}]_{m}}\bm{X}_{e}(\bm{q})\\ \end{array}\right]_{ii}$$

$$\begin{split}\frac{\partial h(\bm{q},\bm{\lambda}^{pre})}{\partial[\bm{q}]_{m}}{}={}&Re\left\{\begin{array}[]{c}\sum\limits_{i=1}^{n_{x}}\overline{([eig(\bm{A}(\bm{q}))]_{i}-[\bm{\lambda}^{pre}]_{i})}\\ \end{array}\right.\\ &\left.\begin{array}[]{c}\times\frac{\partial[eig(\bm{A}(\bm{q}))]_{i}}{\partial[\bm{q}]_{m}}\\ \end{array}\right\}\mbox{ and }\\ [\bm{H}]_{m\ell}{}={}&Re\left\{\begin{array}[]{c}\sum_{i=1}^{n_{x}}\overline{(\frac{\partial[eig(\bm{A}(\bm{q}))]_{i}}{\partial[\bm{q}]_{m}})}(\frac{\partial[eig(\bm{A}(\bm{q}))]_{i}}{\partial[\bm{q}]_{\ell}})\\ \end{array}\right\}\end{split}$$

$$\begin{array}{rl}{\mathbf{p}}_k={\arg }_{\mathbf{p}}\underset{\mathbf{p}}{\min }& m_k(\mathbf{p})\\ \text{mbox}subjectto& ||\mathbf{p}|{|}_2\le {\Delta }_k\end{array}$$

$$m_k(\mathbf{p})=h({\mathbf{q}}_k,{\mathbf{\lambda }}^{pre})+\nabla h({\mathbf{q}}_k,{\mathbf{\lambda }}^{pre}{)}^T\mathbf{p}+{\mathbf{p}}^T{\mathbf{H}}_k\mathbf{p}$$

$${\varphi }_k=\frac{h({\mathbf{q}}_k,{\mathbf{\lambda }}^{pre})-h({\mathbf{q}}_k+{\mathbf{p}}_k,{\mathbf{\lambda }}^{pre})}{m_k(0)-m_k({\mathbf{p}}_k)}.$$

$${\Delta }_{k+1}:={\Delta }_k/4$$

$${\Delta }_{k+1}:=\min \{2{\Delta }_k,\hat{\Delta }\}$$

Full text

I Introduction

Bilinear matrix inequality (BMI) problems frequently arise in system and controller designs [1], e.g., low-authority controller (LAC) designs [2, 3], static output feedback designs for spectral abscissa optimization/$H_{2}$ optimization/$H_{\infty}$ optimization [4, 5, 6], affine fuzzy system designs [7, 8], and observer-based robust controller designs [9]. The advantages of using BMI formulations can be observed in various scenarios. For instance, BMI formulations can avoid a nonsmooth objective function that is hard to handle when spectral abscissa optimization is considered [10]; they may outperform linear matrix inequality (LMI) approaches that can fail to predict the stability of Takagi–Sugeno fuzzy systems [11]; and they can yield less conservative designs than using LMI formulations [12].

While BMI problems are NP-hard [13, 14], BMI solution methods are continuously investigated in the literature because of the advantages derived from using BMI formulations. In [2, 15, 16], path-following methods were proposed in which controller gains were iteratively perturbed to achieve desired performance specifications. The methods were based on the assumption that closed- and open-loop systems were slightly different, i.e., LAC designs were considered. In [4], convex–concave decomposition and linearization methods (CCDM) were combined to address static output feedback problems. After decomposition and linearization, BMI constraints were addressed by solving a sequence of convex semi-definite programming problems. In [5], an inner convex approximation method (ICAM) was proposed as a generalized version of the CCDM. Nonlinear semi-definite programming was considered and a regularization technique was employed to ensure a strict descent search direction. In [17], a Newton-like search method closely related to alternative projection methods was proposed to improve convergence properties.

Alternating minimization (AM) is another popular solution method and has been widely used because of their simplicity and effectiveness [11, 18, 19]. For the AM methods, decision variables are divided into two groups. By fixing one group of variables, the other group of variables forms an LMI problem (LMIP), which is convex and can be solved efficiently. Decision variables in separate groups are then determined alternately during the solving process of LMIPs. Variant versions include iterative LMI (ILMI) methods [7, 20, 21] and the two-step procedure [9]. A few Matlab toolboxes for BMI problems are also available online. For example, LMIRank can be used to solve rank constrained LMI problems [22]. HIFOO employs quasi-Newton updating and gradient sampling to search for solutions. It mainly focuses on fixed-order stabilization and performance optimization problems [23, 24, 25]. PENBMI, commercial software, aims at solving BMI constrained optimization problems or optimization problems that have quadratic cost functions [26, 27].

The aforementioned methods and software packages serve as local optimization approaches to BMI problems. Because a BMI problem is nonconvex, local optima exist and, hence, local optimization approaches may not be able to achieve global optimality. To avoid attaining local optimality, we consider global optimization approaches that employ heuristic algorithms. In [28, 29, 30, 31, 32, 33], branch-and-bound (BB) type methods were proposed. The BB type methods replace bilinear terms with bounded new variables so that a BMI problem can be relaxed into an LMIP. Although being possible to achieve the global optimum, BB type methods can bear a computational burden because the size of the LMIs that must be solved for the lower bound can increase exponentially upon increasing the number of decision variables [21]. In [34], another global optimization approach using generalized Benders decompositions was proposed for BMI problems, but its performance was not evaluated through a number of test problems.

In general, existing BMI solution methods can suffer from at least one of the following five drawbacks or limitations. First, decision variables are expressed solely in a vector form, e.g., some BB type methods. By contrast, a matrix form is more convenient in control problems [35]. Second, solution methods are originally designed to fit particular problem structures. In some situations, applying developed methods to other problem structures, if not impossible, requires extra efforts to reformulate the problem, e.g., some AMs and ILMI methods. In other situations, solution methods cannot be applied to problems that do not have the intended structures, e.g., path-following methods. Third, prior derivations such as approximations or decompositions must be performed before algorithms are applied, e.g., the CCDM and ICAM, and these derivations can be cumbersome and sometimes heuristic. Fourth, only local optimization is performed while BMI problems inherently have multiple local optima. Finally, to the best of our knowledge, existing BMI solution methods cannot address multiobjective optimization problems (MOPs) in which a set of Pareto optimal solutions is of interest rather than the global optimal solution.111The ability to solve MOPs constrained by BMIs is worth further investigation because MOPs naturally and frequently arise in engineering problems [36, 37, 38, 39]. Solving an MOP, yielding an approximate Pareto front (APF) and Pareto optimal set, can provide a system designer with a broad perspective on optimality. The resulting APF can clearly illustrate how one objective affects the others, and the obtained Pareto set allows the designer to make a posterior decision, i.e., selecting design parameters after a set of promising candidates is available [40, 41]. In general, a posterior decision is preferred to a prior decision because more information has been used before the decision making [42].

To avoid the aforementioned five drawbacks or limitations, we propose a method of reduction of variables (MRV). The method consists of a principle of variable classification, a transformation of the BMI problem, and a hybrid multiobjective immune algorithm (HMOIA) that solves the problem derived from the transformation. Internal and external variables are coined and used to denote all the decision variables involved. The internal variable can represent a set of matrix variables, which is convenient in controller designs. To develop a general-purpose solution method, we consider possible multiple objectives in BMI problems and assume no particular problem structures. This yields a framework that addresses feasibility problems, single-objective optimization problems (SOPs), and MOPs constrained by BMIs in a unified manner. The developed HMOIA is a hybrid because it employs stochastic and deterministic mechanisms to determine the external and internal decision variables, respectively. The stochastic mechanism allows for global exploration of the entire solution space. By applying the HMOIA to BMI problems, few prior derivations, involving only variable classification and simple problem transformation, are required. Limited derivations render the proposed method suitable for various BMI problems.

To verify the effectiveness of the MRV, we used a series of test problems in our simulations [6, 11, 7, 8, 9, 29, 2, 15]. For feasibility problems, while different solution methods were developed to address various BMI problems, the MRV was able to find a solution with $100\%$ success rates in a unified manner. In spectral abscissa optimization, the MRV outperformed existing methods in $73\%$ of selected benchmark problems in terms of the minimum value or mean value. The MRV achieved better levels of performance than existing methods in $27.5\%$ and $47.8\%$ of selected $H_{2}$ and $H_{\infty}$ optimization problems, respectively, while it yielded similar performance in the remaining problems. As shown in [4] and [5], the CCDM and ICAM were relatively robust compared with other existing solution methods. We illustrated that the MRV was able to find solutions to certain problems in which these two robust methods failed or made little progress towards a local solution.

The main contributions of this paper are as follows. We propose a novel global optimization approach to BMI problems, which has not been fully investigated compared to local optimization approaches. This approach can combat a few drawbacks existing BMI solution methods can suffer from: using inconvenient variable expression, being confined to particular problem structures, requiring heuristic or cumbersome prior derivations, or being incapable of addressing multiple objectives. When the proposed classification principle is applicable, we provide a unified formulation that facilitates generating solutions to feasibility problems, SOPs, and MOPs constrained by BMIs. To the best of our knowledge, this is the first study that provides such a unified framework. We perform related analysis and validate the proposed MRV through a large number of benchmark problems, showing that the proposed methodology can outperform existing solution methods in many of these BMI problems.

The rest of this paper is organized as follows. Section II describes the problem formulation and the principle of variable classification. In Section III, preliminaries to our algorithm development are examined, including analysis of problem transformations. Section IV presents the HMOIA and hence, the MRV. Simulation results are given in Section V. Finally, Section VI concludes this paper.

II Problem Formulation and Variable Classification

In this section, we investigate system and control designs that are formulated as BMI problems, and propose a classification principle for decision variables that facilitates solution search. Under our framework, the associated cost function can be a vector-valued function, a scalar function, or a constant, depending on the number of objectives involved. By using the classification principle, decision variables in BMIs are classified into two types, the internal and external variables. Design examples are presented to illustrate how to use the proposed principle of variable classification.

The following notation and terminology are used throughout this study. Let $\mathbb{R}$ and $\mathbb{C}$ be the sets of real and complex numbers, respectively. For a scalar $b\in\mathbb{C}$, $\overline{b}$ denotes the complex conjugate of $b$. Let $[\bm{a}]_{i}$ and $[\bm{A}]_{ij}$ denote the $i$th entry of the vector $\bm{a}$ and the $(i,j)$th entry of the matrix $\bm{A}$, respectively. For two vectors $\bm{a}$ and $\bm{b}$, $\bm{a}\leq\bm{b}$ is interpreted as $[\bm{a}]_{i}\leq[\bm{b}]_{i}$ for all $i$. If $\bm{P}>0$, then $\bm{P}$ is symmetric and positive-definite. Similarly, $\bm{P}<0$ implies that $\bm{P}$ is symmetric and negative-definite. For a square matrix $\bm{A}$, $eig(\bm{A})$ represents the vector of all eigenvalues of $\bm{A}$ placed in a prescribed manner, and $eig\{\bm{A}\}$ represents the set of all eigenvalues of $\bm{A}$. The mark “$\star$” is used to denote the induced symmetry, e.g., $(\bm{PA},\star)=\bm{PA}+\bm{A}^{T}\bm{P}^{T}$ and

[TABLE]

If $\bm{f}:\Omega\rightarrow\mathbb{R}^{N}$ is a vector-valued function, then the MOP

[TABLE]

is interpreted as vector optimization in which Pareto optimality is adopted. The domain $\Omega$ lies in the Euclidean space $\mathbb{R}^{M}$ for some positive integer $M$. The associated terminology is presented as follows [43, 44, 45].

Definition 1** (Pareto dominance)**

In the decision variable space of  (1), a point $\bm{\omega}^{\prime}\in\Omega$ dominates another point $\bm{\omega}^{\prime\prime}\in\Omega$ if the conditions $[\bm{f}(\bm{\omega}^{\prime})]_{i}\leq[\bm{f}(\bm{\omega}^{\prime\prime})]_{i},i=1,2,...,N,$ hold true and at least one inequality is strict. In this case, we denote $\bm{\omega}^{\prime}\preceq_{\bm{f}}\bm{\omega}^{\prime\prime}$ and $\bm{f}(\bm{\omega}^{\prime})\preceq\bm{f}(\bm{\omega}^{\prime\prime})$. A point that is not dominated by other points is termed a nondominated point.

Definition 2** (Pareto optimal set)**

The Pareto optimal set $\mathcal{P}^{*}$ of (1) is defined as the set of all nondominated points, i.e.,

[TABLE]

Definition 3** (Pareto front)**

The Pareto front (PF) of (1) is defined as the image of the Pareto optimal set through the mapping $\bm{f}$, i.e., $\bm{f}(\mathcal{P}^{*})$ represents the PF.

In our BMI-based design problems, we use

[TABLE]

to represent a BMI, where $\mathcal{BMI}(\cdot)$ is a matrix function, and $\bm{\alpha}$ and $\bm{X}$ are the variables. The inequality $\mathcal{BMI}(\bm{\alpha},\bm{X})<0$ becomes an LMI in the variable $\bm{\alpha}$ given $\bm{X}$ or in the variable $\bm{X}$ given $\bm{\alpha}$. If more than one BMI are involved, then the notation $\mathcal{BMI}(\bm{\alpha},\bm{X})$ represents a block-diagonal matrix such that $\mathcal{BMI}(\bm{\alpha},\bm{X})<0$ consists of all the BMIs.

To consider optimal designs in a unified framework, we add an objective function $\bm{\mathcal{F}}(\cdot)$ to (2). From the perspective of algebra, there is no difference between $\bm{\alpha}$ and $\bm{X}$ in (2) because they are just two coupled variables in the BMI $\mathcal{BMI}(\bm{\alpha},\bm{X})<0$. However, to create a solution method, we assume $\bm{\mathcal{F}}$ is a function of $\bm{\alpha}$. The resulting BMI-based MOP can be expressed as

[TABLE]

where $\bm{\alpha}$ is distinguished from $\bm{X}$ by using the following classification principle.

Principle of Variable Classification:

Upper and lower bounds on the entries of variables in $\bm{\alpha}$ are available or can be obtained. Square matrix variables in $\bm{\alpha}$, if any, do not have constraints on definiteness, i.e., positive or negative definiteness. 2. 2.

Bounds on entries of variables in $\bm{X}$ are unavailable. 3. 3.

The objective function $\bm{\mathcal{F}}$ can be expressed solely in terms of $\bm{\alpha}$. 4. 4.

The size of $\bm{\alpha}$ should be as small as possible.

When the classification principle is applicable, we term the variables $\bm{\alpha}$ and $\bm{X}$ the external and internal decision variables, respectively. In our principle, $\bm{\alpha}$ represents those variables (scalar and/or matrix variables) in a BMI problem that have bounds on entries. These bounds are mostly inherent from physical constraints or can be readily assigned mathematically. The remaining variables (scalar and/or matrix variables) are included in $\bm{X}$. They generally do not have upper and lower bounds on their entries, but there can be constraints related to positive or negative definiteness imposed on matrix variables in $\bm{X}$. The definiteness associated with matrix variables in $\bm{X}$ is required to ensure the system stability, which mainly distinguishes $\bm{X}$ from $\bm{\alpha}$. A typical external variable can include controller gains and/or system parameters. By contrast, matrix variables related to the Lyapunov theory are classified as the internal variable because bounds on the entries of these matrix variables are unavailable in practice.

The condition in which $\bm{\mathcal{F}}$ is not a function of $\bm{X}$ does not yield a restricted problem formulation. For instance, if $[\bm{\mathcal{F}}]_{i}=g(\bm{X})$ is encountered, we may introduce a slack variable $\eta$, impose the constraint $g(\bm{X})\leq\eta$, and assign $[\bm{\mathcal{F}}]_{i}:=\eta$. In this way, the objective function becomes the one with $\bm{\alpha}$ as the only variable. Finally, it will be shown that the BMI-constrained problem in (3) can be reduced to an unconstrained problem in which $\bm{\alpha}$ is the only decision variable. Therefore, a smaller size of $\bm{\alpha}$ means the fewer number of decision variables in the unconstrained problem, which explains why we keep the size of $\bm{\alpha}$ as small as possible in the classification principle.

When $\bm{\mathcal{F}}(\cdot)$ is a constant function, it is understood that the BMI problem in (3) is interpreted as a feasibility problem. Otherwise, an SOP (or MOP) is considered if $\bm{\mathcal{F}}(\cdot)$ is a scaler-valued (or vector-valued) function. For a feasibility problem, it is desired to determine whether or not there exists a point $(\bm{\alpha},\bm{X})$ satisfying the matrix inequality $\mathcal{BMI}(\bm{\alpha},\bm{X})<0$. If such a point exists, then the problem is feasible and any point that satisfies the matrix inequality is a solution (or a feasible point). For an SOP, it is desired to search for a feasible point that achieves the minimum value of the objective function. When an MOP is considered, the associated optimality is interpreted as Pareto optimality. In that case, the Pareto optimal set is to be determined.

To illustrate how to use the principle of variable classification, we examine a few design examples as follows.

II-A Feasibility Problems

Stability Test (ST): Consider a T–S fuzzy system [46]

[TABLE]

It can be shown that the system in (4) is stable if there exist $\tau_{\ell ij}\geq 0$ and $\bm{P}_{i}>0$ such that [11]

[TABLE]

are satisfied. According to the classification principle, the external variable cannot include matrix variables that have a constraint on definiteness. Because $\bm{P}_{i}>0,i=1,2,$ are positive-definite matrix variables, they must be included in the internal variable $\bm{X}$; to yield a BMI problem, the remaining variables $\tau_{\ell ij}$ are included in the external variable $\bm{\alpha}$. The feasibility problem in (5) can then be expressed as $\mathcal{BMI}(\bm{\alpha},\bm{X})<0$ in which $\bm{\alpha}=(\tau_{112},\tau_{121},\tau_{212},\tau_{221})$ and $\bm{X}=(\bm{P}_{1},\bm{P}_{2})$.

II-B Single-objective Optimization Problems

Linear Parameter-varying Systems (LPVS): Consider a linear time-varying system [47, 48]

[TABLE]

where

[TABLE]

The $\varsigma$ represents a design parameter. The system in (6) is stable if there exist $\delta_{i}$ and $\bm{P}_{i}$ satisfying [29, 49]

[TABLE]

For a fixed $\varsigma$, (7) is a BMI in the variables $(\delta_{1},\delta_{2})$ and $(\bm{P}_{1},\bm{P}_{2})$. To find the largest value of $\varsigma$ yielding a stable system, we can solve

[TABLE]

Based on the principle of variable classification, $\bm{P}_{1}$ and $\bm{P}_{2}$ are positive-definite and must be included in the internal variable $\bm{X}$; to have $\mathcal{BMI}(\bm{\alpha},\bm{X})<0$ as an LMI problem for a fixed $\bm{\alpha}$, we are forced to include all the remaining variables in the external variable. We thus have $\bm{\alpha}=(\varsigma,\delta_{1},\delta_{2})$, $\bm{X}=(\bm{P}_{1},\bm{P}_{2})$, and $\bm{\mathcal{F}}(\bm{\alpha})=-\varsigma$. The negative sign in $\bm{\mathcal{F}}$ has been added for the conversion of (8) to the minimization form of (3).

II-C Multiobjective Optimization Problems

For a sparse linear constant output-feedback design, the BMI problem [2, 4]

[TABLE]

can be formulated, where $\bm{A}_{\bm{F}}=\bm{A}+\bm{BFC}$, $\sigma>0$ represents a prescribed weighting coefficient, and $\beta$ represents the decay rate. The SOP in (9) is interpreted as determining the controller gain $\bm{F}$ so that the decay rate $\beta$ is maximized and $\bm{F}$ is kept as much sparse as possible. One drawback of considering the single-objective formulation is that there is no rule that can be used to assign the value of $\sigma$, which affects the values of $\beta$ and $\bm{F}$. In practice, a system designer selects an arbitrary value of $\sigma$ and accepts the resulting gain $\bm{F}$. To avoid such heuristic assignment for $\sigma$, we can consider a multiobjective formulation that addresses two objectives in separate dimensions [42, 50]:

[TABLE]

According to the classification principle, $\bm{P}$ is positive-definite and hence, included in the internal variable $\bm{X}$; to yield a BMI problem, the remaining variables must be included in the external variable $\bm{\alpha}$. Referring to (3), we have $\bm{\alpha}=(\beta,\bm{F})$, $\bm{X}=\bm{P}$, and $\bm{\mathcal{F}}(\bm{\alpha})=[-\beta\;\sum_{i}\sum_{j}|[\bm{F}]_{ij}|]^{T}$. Once (10) has been solved, an approximate Pareto front (APF) can be obtained and the system designer can select an appropriate $\bm{F}$ based on the information provided by the APF.

The proposed classification principle is based on the basic properties of BMIs represented by $\mathcal{BMI}(\bm{\alpha},\bm{X})<0$. Because BMIs are nonlinear and have possibly several local optima when optimization is involved, any deterministic algorithms can be trapped locally. To remedy this problem, stochastic algorithms can be used. Because $\mathcal{BMI}(\bm{\alpha},\bm{X})<0$ is a BMI, $\mathcal{BMI}(\bm{\alpha},\bm{X})<0$ becomes an LMI in the variable $\bm{X}$ for a fixed value of $\bm{\alpha}$. For LMIs, it is well-known that deterministic algorithms such as interior-point methods are suitable for solving them efficiently. These arguments suggest that variables in BMI problems be classified into two groups so that a hybrid algorithm combining stochastic and deterministic search engines can be applied.

To integrate stochastic and deterministic search schemes, we first explore the variable space of $\bm{\alpha}$ (external exploration) so that an LMIP $\mathcal{BMI}(\bm{\alpha},\bm{X})<0$ in the variable $\bm{X}$ can be obtained. Once $\bm{\alpha}$ is determined, the associated variable space of $\bm{X}$ can then be searched internally and efficiently because of the convexity. This explains why $\bm{\alpha}$ and $\bm{X}$ are termed external and internal variables, respectively. Since $\bm{X}$ is relevant to the feasibility but irrelevant to the objective values, this internal variable can be considered hidden from the external search if information about the feasibility is extracted properly. Therefore, we may reduce the original problem with variables $\bm{\alpha}$ and $\bm{X}$ to a simpler problem with only the variable $\bm{\alpha}$, and then transform the resulting problem into another form that is convenient for addressing the feasibility condition.

III Preliminaries to Algorithm Development

This section discusses the reduction and transformation, and other preliminary results that are helpful in later development of the hybrid algorithm. The section is divided into three subsections: Section III-A focuses on theorems that transform the BMI problem in (3) into an unconstrained problem with fewer decision variables; Section III-B presents a solution method related to pole placement problems; and Section III-C describes an algorithm that reduces the population density of the HMOIA.

III-A Reduction and Equivalence Theorems

Theorems in this subsection lead to an optimization problem that has a simpler form than (3). By using the theorems, the number of decision variables in (3) can be reduced, and the associated problem can be further transformed into an unconstrained optimization problem. Although we adopt multiobjective formulations in the following discussions, the established results remain true when an SOP or a feasibility problem is considered.

Consider the eigenvalue problem (EVP)

[TABLE]

Because the constraint $\mathcal{BMI}(\bm{\alpha},\bm{X})<0$ is a BMI, the EVP in (11) is convex in the variables $\lambda$ and $\bm{X}$ given the value of $\bm{\alpha}$. (For a fixed value of $\lambda$, the EVP can thus be solved by interior-point methods.) In (11), we denote $(\lambda^{*}(\bm{\alpha}),\bm{X}^{*}(\bm{\alpha}))$ as the pair that achieves the minimum. Both $\lambda^{*}(\bm{\alpha})$ and $\bm{X}^{*}(\bm{\alpha})$ are regarded as a function of $\bm{\alpha}$. The following lemma relates the value of $\lambda^{*}(\bm{\alpha})$ to the feasibility of (3).

Lemma 1

The BMI problem in (3) is feasible if and only if an $\tilde{\bm{\alpha}}$ exists such that the value of $\lambda^{*}(\tilde{\bm{\alpha}})$ in (11) is negative, i.e., $\lambda^{*}(\tilde{\bm{\alpha}})<0$.

Proof:

It can be readily verified by a slight modification of the proof in Lemma 1 of [43] or [51]. ∎

The following theorem follows from using Lemma 1.

Theorem 1** (Reduction Theorem)**

There exists a pair $(\tilde{\bm{\alpha}},\tilde{\bm{X}})$ that is Pareto optimal in (3) if and only if (denoted by $\Leftrightarrow$) $\tilde{\bm{\alpha}}$ is Pareto optimal in

[TABLE]

Proof:

We first prove necessity ($\Rightarrow$). By Lemma 1, we have $\lambda^{*}(\tilde{\bm{\alpha}})<0$ and hence, $\tilde{\bm{\alpha}}$ is a feasible point of (12). Let us proceed by contraposition. Suppose that there exists an $\bm{\alpha}^{\prime}$ dominating $\tilde{\bm{\alpha}}$ in (12), i.e.,

[TABLE]

However, the conditions in (13) implies that

[TABLE]

which yields a contradiction.

To prove sufficiency ($\Leftarrow$), we again use contraposition. Suppose that there exists a pair $(\bm{\alpha}^{\prime},\bm{X}^{\prime})$ dominating $(\tilde{\bm{\alpha}},\tilde{\bm{X}})$, i.e.,

[TABLE]

By Lemma 1, the conditions in (14) are equivalent to those in (13), which implies that $\tilde{\bm{\alpha}}$ is not Pareto optimal in (12). However, this contradicts the Pareto optimality of $\tilde{\bm{\alpha}}$. ∎

According to Theorem 1, the BMI problem in (3) with $\bm{\alpha}$ and $\bm{X}$ as the decision variables can reduce to (12) with $\bm{\alpha}$ as the only decision variable. That is, the number of decision variables is reduced, which explains why Theorem 1 is termed the Reduction Theorem. In the theorem, the expression “Pareto optimal” is replaced by “feasible” if we consider a feasibility problem. In this case, the theorem is exactly the same as Lemma 1. Similarly, we replace “Pareto optimality” with conventional optimality when an SOP is encountered.

By Theorem 1, we can solve (12) for a BMI-based design in place of (3). We further consider an unconstrained problem that is equivalent to (12).

Theorem 2** (Equivalence Theorem)**

Let

[TABLE]

where $\max\{0,\lambda^{*}(\bm{\alpha})\}$ represents the maximum element in the set $\{0,\lambda^{*}(\bm{\alpha})\}$. A point $\tilde{\bm{\alpha}}$ is Pareto optimal in (12) if and only if (denoted by $\Leftrightarrow$) $\tilde{\bm{\alpha}}$ satisfies the condition $\max\{0,\lambda^{*}(\tilde{\bm{\alpha}})\}=0$ and is Pareto optimal in

[TABLE]

Proof:

We prove necessity ($\Rightarrow$). Since $\tilde{\bm{\alpha}}$ is Pareto optimal in (12), we have $\lambda^{*}(\bm{\alpha})<0$ and thus $\tilde{\bm{\alpha}}$ satisfies the condition $\max\{0,\lambda^{*}(\tilde{\bm{\alpha}})\}=0$. We use contraposition. Suppose that $\tilde{\bm{\alpha}}$ is not Pareto optimal in (16). There must exist an $\bm{\alpha}^{\prime}$ such that $\bm{\alpha}^{\prime}\preceq_{\tilde{\bm{\mathcal{F}}}}\tilde{\bm{\alpha}}$, yielding $\max\{0,\lambda^{*}(\bm{\alpha}^{\prime})\}=0$. However, this implies that the conditions in (13) hold true, i.e., $\tilde{\bm{\alpha}}$ is not Pareto optimal in (12), which yields a contradiction.

We now prove sufficiency ($\Leftarrow$) and again use contraposition. Suppose that an $\bm{\alpha}^{\prime}$ exists such that the conditions in (13) hold true. This implies $\bm{\alpha}^{\prime}\preceq_{\tilde{\bm{\mathcal{F}}}}\tilde{\bm{\alpha}}$, which yields a contradiction. ∎

Theorem 2 is termed the Equivalence Theorem because it establishes an equivalence relation between (16) and (12). According to Reduction and Equivalence Theorems, we can solve the unconstrained problem in (16) that has fewer decision variables than the original BMI problem in (3).

III-B Levenberg–Marquardt Method

Pole placement problems occur frequently in controller designs [52, 53]. In this subsection, we investigate a trust region Levenberg–Marquardt method that can be used for pole placement when the system matrix $\bm{A}+\bm{BFC}$ is encountered. In this situation, matrices $\bm{A}$, $\bm{B}$, and $\bm{C}$ are known, and $\bm{F}$ is the design parameter that must be determined.

Suppose that $\bm{A}\in\mathbb{R}^{n_{x}\times n_{x}}$ and $\bm{F}\in\mathbb{R}^{n_{u}\times n_{y}}$ in which $n_{x},n_{y},$ and $n_{u}$ represent the dimensions of the state vector, physical output, and control input, respectively. To facilitate the following discussions, we reshape the gain matrix $\bm{F}$ into an $n_{u}n_{y}\times 1$ vector $\bm{q}$, and denote $\bm{A}(\bm{q})=\bm{A}+\bm{BFC}\in\mathbb{R}^{n_{x}\times n_{x}}$ and

[TABLE]

where $eig(\bm{A}(\bm{q}))$ represents the vector of eigenvalues of $\bm{A}(\bm{q})$ and $\bm{\lambda}^{pre}$ is a prescribed vector of poles. The entries of the vector $eig(\bm{A}(\bm{q}))$ in (17) is placed in a way that the minimum norm is achieved. The associated pole placement problem can be formulated as

[TABLE]

which is an unconstrained nonlinear least squares problem. To apply the trust region Levenberg–Marquardt method to solve (18), we need the first partial derivatives and an approximate Hessian matrix of $h(\bm{q},\bm{\lambda}^{pre})$ in (17). For the vector of eigenvalues $eig(\bm{A}(\bm{q}))$, let $\bm{X}_{e}(\bm{q})\in\mathbb{C}^{n_{x}\times n_{x}}$ be the matrix consisting of the associated eigenvectors such that

[TABLE]

We have [54, 55]

[TABLE]

for $i=1,2,...,n_{x}$, and $m=1,2,...,n_{u}n_{y}$. The first partial derivatives $\partial h(\bm{q},\bm{\lambda}^{pre})/\partial[\bm{q}]_{m}$ and approximate Hessian matrix $\bm{H}$ can be expressed as

[TABLE]

respectively. The Levenberg–Marquardt algorithm for the pole placement problem in (18) is described as follows.222The reader can refer to Theorems 1–4 in [56] or Theorems 4.8, 4.9, and 6.4 in [57] for the convergence analysis of the trust region method.

Trust Region Levenberg–Marquardt Algorithm [56, 57]

 

Given $\hat{\Delta}>0$, $\Delta_{0}\in(0,\hat{\Delta})$, and $\eta\in[0,1/4)$

For $k=1,2,...$

• Evaluate $\bm{p}_{k}$ by solving

[TABLE]

where $\Delta_{k}$ represents the current trust region radius, and

[TABLE]

with entries of $\nabla h(\bm{q}_{k},\bm{\lambda}^{pre})$ and $\bm{H}_{k}$ defined in (19).

• Evaluate

[TABLE]

• If $\phi_{k}<1/4$

[TABLE]

Else

• If $\phi_{k}>3/4$ and $||\bm{p}_{k}||_{2}=\Delta_{k}$

[TABLE]

Else

[TABLE]

End If

End If

• If $\phi_{k}>\eta$

[TABLE]

Else

[TABLE]

End If

End For

III-C Density Reduction Algorithm

When an evolutionary algorithm searches for Pareto optimal solutions to an MOP, less crowded points must be preserved so that population diversity can be ensured. To this end, we estimate the density of current population and remove points that lie in a dense region. We denote $\mathcal{A}(t_{c})$ as the current population with the cardinality $|\mathcal{A}(t_{c})|$, $N_{nom}$ as the nominal size of the population, and $\bm{\alpha}$ as an element of $\mathcal{A}(t_{c})$. Suppose that $\tilde{\bm{\mathcal{F}}}(\bm{\alpha})\in\mathbb{R}^{N+1}$. The process of removing points from a dense region is termed density reduction, which can be performed by the following algorithm modified from [58].

**Density Reduction Algorithm

 

While** $|\mathcal{A}(t_{c})|>N_{nom}$ **do

**Evaluate

[TABLE]

for all $\bm{\alpha}_{i}\in\mathcal{A}(t_{c})$ and $j=1,2,...,N$, where

[TABLE]

Evaluate

[TABLE]

for all $\bm{\alpha}_{i}\in\mathcal{A}(t_{c})$.

Remove the element $\bm{\alpha}$ that yields the least $(\bm{\alpha})_{av}$ from $\mathcal{A}(t_{c})$ and thus the size of $\mathcal{A}(t_{c})$ is reduced by one.

**End While

 **

IV Proposed Algorithm

This section presents the HMOIA used to solve (16). The algorithm is a hybrid because it integrates both stochastic and deterministic search schemes. For example, the Levenberg–Marquardt algorithm, density reduction algorithm, and interior-point methods are deterministic algorithms, while artificial immune systems used as the underlying structure of the HMOIA are stochastic search methods. There are a few reasons why the immune search scheme was adopted in our main algorithm structure. First, its potential to provide novel solutions has been illustrated in several studies [59, 60, 61]. Second, the immune search scheme is robust and outperforms some existing MOEAs or at least performs equally well in most benchmark MOPs [58]. (In [58], the MO immune algorithm was compared to PAES, PESA, NSGA-II [62], SPEA2, MOEA/D [63], and ACSAMO in terms of convergence, diversity, uniformity, and coverage.) Finally and most importantly, the artificial immune system is a “highly parallel intelligent system” [64, 65, 42] and thus a parallel computation scheme can be readily developed. This is useful for solving BMI constraints that requires a large amount of computational power in general. Despite these reasons, it is worth mentioning that other advanced MOEAs can also be adopted if modified properly.333Proper modifications may include incorporation of pole-placement techniques into the search engine and a design of a mechanism that ensures legitimate pole placement.

The pseudocode of the HMOIA is presented as follows.

**Pseudocode of the Proposed HMOIA

 

Input: **MOP in (16)

Prescribe bounds on the external variable and initialize the population

Evaluate the objective function

Remove dominated points

Let $t_{c}:=1$

While $t_{c}\leq t_{max}$ do

• If $N=0$ and $\exists\bm{\alpha}\in\mathcal{A}(t_{c})$ such that $\tilde{{\bm{\mathcal{F}}}}(\bm{\alpha})=0$

Let $t_{c}:=t_{max}$

Else

• Perform the hyper-mutation operation
• Evaluate the objective function
• Update population

End If

• Let $t_{c}:=t_{c}+1$

**End While

**Remove points $\bm{\alpha}$ that have $[\tilde{\bm{\mathcal{F}}}(\bm{\alpha})]_{N+1}>0$

Remove dominated points

**Output: **Approximate Pareto optimal solutions and Pareto front

In the following subsections, we elaborate key steps of the algorithm and summarize the MRV.

IV-A Prescribe Bounds on External Variable and Initialize the Population

To specify the range of interests, we prescribe bounds for the external variable $\bm{\alpha}$ in (16). Entries of $\bm{\alpha}$ can be generated pointwisely over prescribed bounds or recovered collectively from a given vector of eigenvalues $\bm{\lambda}^{pre}$ described in (17) and (18). For example, if a range $[\alpha^{min}_{i},\alpha^{max}_{i}]$ is prescribed, then the $i$th entry of $\bm{\alpha}$ can be generated uniformly at random over the range. Otherwise, if a range $[-\sigma^{min},0]\times[-\omega^{max},\omega^{max}]$ is given, we can randomly generate

[TABLE]

where complex entries of $\bm{\lambda}^{pre}$ occur in conjugate pairs, and then recover the entries of $\bm{\alpha}$ from $\bm{q}^{*}(\bm{\lambda}^{pre})$ defined in (18) using the trust region Levenberg–Marquardt algorithm presented in Section III-B.

After specifying the range, we initialize the population: assign the nominal population size $N_{nom}$ and the maximum population size $N_{max}$, and generate initial population

[TABLE]

The basic structure of artificial immune algorithms in [45] and [58] is adopted. During the evolutionary process that mimics operations in immune systems, the population size changes over the iteration but remains below $N_{max}$, nondominated points are maintained, and dominated points are removed from the population.

In general, large bounds should be assigned to provide a spacious search space, but such a spacious space can yield ineffective search when Pareto optimal solutions have most entries that are close to zero. To manage possible ineffectiveness, we divide the search spaces $[\alpha^{min}_{i},\alpha^{max}_{i}]$ or $[-\sigma^{min},0]\times[-\omega^{max},\omega^{max}]$ into several subspaces

[TABLE]

where $\kappa_{s}\in(0,1]$ with $s=1,2,...,S$. Population initialization is thus modified accordingly. We either generate $[\bm{\alpha}]_{i}\in[\kappa_{s}\alpha^{min}_{i},\kappa_{s}\alpha^{max}_{i}]$ pointwisely or recover entries of $\bm{\alpha}$ collectively from

[TABLE]

IV-B Evaluate the Objective Function

Function evaluation for the objective function $\tilde{\bm{\mathcal{F}}}$ defined in (15) can be divided into two parts. The first part addresses the evaluation of $\bm{\mathcal{F}}(\bm{\alpha})$. If $\bm{\mathcal{F}}$ is explicitly expressed as a function of $\bm{\alpha}$, then the evaluation is simply the substitution of $\bm{\alpha}$ into $\bm{\mathcal{F}}$; otherwise, deterministic algorithms are employed to evaluate $\bm{\mathcal{F}}(\bm{\alpha})$. For example, if $\bm{\alpha}$ represents a controller gain of a linear control system and $\bm{\mathcal{F}}(\bm{\alpha})$ denotes the associated $H_{\infty}$ norm, then $\bm{\mathcal{F}}(\bm{\alpha})$ must be evaluated using deterministic algorithms. The second part addresses the evaluation of $\lambda^{*}(\bm{\alpha})$. Because this evaluation is related to solving LMIs, deterministic algorithms such as interior-point methods can be used.

IV-C Remove Dominated Points

By removing dominated points from the population, nondominated points are maintained. Preserving nondominated points is an important operation that relates to the convergence of the algorithm. Other operations such as the hyper-mutation and population update that guide the population towards the Pareto optimal set are important to the algorithm convergence as well.

IV-D Perform Hyper-mutation Operation

Let $\mathcal{A}(t_{c})$ and $|\mathcal{A}(t_{c})|$ denote the current population and the associated population size, respectively. For two vectors $\bm{a}$ and $\bm{b}$, $\bm{a}\oplus\bm{b}$ denotes a random and pointwise combination of entries of $\bm{a}$ and $\bm{b}$, i.e., $[\bm{a}\oplus\bm{b}]_{i}$ can be either $[\bm{a}]_{i}$ or $[\bm{b}]_{i}$ with equal probability. For a hyper-mutation operation, new points $\bm{\alpha}^{j}_{i}$ are generated by

[TABLE]

for all $\bm{\alpha}_{i}\in\mathcal{A}(t_{c})$, where $rand$ and $L^{j}$ are independent random numbers chosen from $(0,1)$, and entries of $\bm{\alpha}^{j}$ are generated pointwisely or recovered collectively in the same way described at the population initialization. The operation in (28) is performed $R(t_{c})=\llcorner N_{max}/|\mathcal{A}(t_{c})|\lrcorner$ times for each $i$, where $\llcorner\cdot\lrcorner$ represents the floor function. This operation can be interpreted as follows: an $\bm{\alpha}_{i}$ in $\mathcal{A}(t_{c})$ is cloned $R(t_{c})$ times, and then all these cloned points mutate to produce points $\bm{\alpha}^{j}_{i},j=1,2,...,R(t_{c})$.

IV-E Update Population

Updating the population consists of addition and removal operations. After the hyper-mutation operation, $R(t_{c})\times|\mathcal{A}(t_{c})|$ points are newly generated and added to the population. To keep a manageable size of the population, we remove infeasible points, dominated points, or nondominated points in order if necessary. Fig. 1 shows a removal procedure that reduces the size of $\mathcal{A}(t_{c})$ to $N_{nom}$. In the procedure, infeasible points $\bm{\alpha}$ that have $\lambda^{*}(\bm{\alpha})>0$ are gradually removed from the population. After the removal, if the population size is still greater than its nominal size, then dominated points are removed from the population randomly and iteratively. If $|\mathcal{A}(t_{c})|$ is still greater than $N_{nom}$, then we remove nondominated points using the density reduction algorithm described in Section III-C.

IV-F Parameter Selection

The values of parameters $N_{max}$, $N_{nom}$, and $t_{max}$ can affect algorithm performance. In general, larger values of them yield a better level of performance if the complexity is not a concern [58, 66, 67, 68]. This is because larger values of $N_{max}$ and $N_{nom}$ mean that more computational resources are employed to explore the search space in each iteration, and a larger value of $t_{max}$ corresponds to more exploration time. When the values of the parameters exceed certain thresholds, mature convergence is attained and further improvement can be hardly observed. Some practitioners suggest that large values of $N_{max}$, $N_{nom}$, and $t_{max}$ be set first, and then these values be lowered gradually until unacceptable results are obtained. However, this practice suffers from two drawbacks. First, the notion of unacceptable results is vague. Second, computational time can be a cost, and repeating the whole search with different parameters becomes costly.

Although parameters should be set differently in different problems for better performance, we selected the same parameters in our simulations when solving all the BMI problems. The main reason of such selection is that it is difficult to define the “optimal” values for parameters in consideration of the performance, complexity, and computational time. Despite of using the same parameters, the simulation results still provide a proof of concept that the proposed methodology can outperform existing design approaches in most benchmark BMI-based design problems.

IV-G Method of Reduction of Variables

For a system or control design problem constrained by BMIs, the proposed method of solution can be realized by the following steps:

S1) Classify the decision variables into the internal and external variables using the Principle of Variable Classification.

S2) Transform the BMI problem in (3) into its equivalent form in (16).

S3) Apply the proposed HMOIA to solve (16).

After these three steps, the associated system or controller can be constructed based on the obtained solution(s).

Remark 1

There are two circumstances in which the proposed methodology can fail: the classification principle is not applicable; and same eigenvalues are assigned when pole placement is performed. The first situation may occur when physical or mathematical bounds on coupled decision variables cannot be obtained or readily prescribed. Here is an example:

[TABLE]

where $\bm{P}$ and $\bm{F}$ are decision variables, yielding a BMI problem. In (29), we cannot assign $\bm{\alpha}=\bm{P}$ or $\bm{\alpha}=\bm{F}$ because both of them do not have inherent bounds on their entries and have a constraint of positive definiteness. The classification principle is thus not applicable. Fortunately, although the proposed principle is not valid, we rarely encounter this type of problem such as (29) in system and control designs. When a controller design problem is considered, $\bm{F}$ generally relates to a controller gain and physical constraints do not yield a requirement of positive or negative definiteness on $\bm{F}$. In fact, $\bm{F}$ may not even be a square matrix in practice. The second situation is related to the differentiability of eigenvalues. Eigenvalues are differentiable only if they are distinct, and the condition of differentiability is used when pole placement is performed. Since eigenvalues are randomly assigned in our algorithm, there is little chance that two eigenvalues are the same. Even if their values are slightly different, the trusted region algorithm described in Section III-B can still work [56]. Therefore, this situation does not impose a serious restriction on the applicability of our methodology either.

Remark 2

In our framework, SOPs and MOPs are addressed in a unified manner. For an illustrative purpose, we examine the problem

[TABLE]

which can represent an SOP or MOP depending on the dimension of $f(\cdot)$. The problem in (30) can then be transformed into the MOP

[TABLE]

If an SOP is considered in (30), then the resulting problem in (31) is a 2-D MOP. If an MOP with two objectives is considered in (30), then the resulting problem in (31) becomes a 3-D MOP. For either case, (31) is regarded as an MOP and can be solved by our hybrid algorithm, producing a solution set. If any points in the obtained solution set yield a nonzero value of the final objective, i.e., $\max\{g(\bm{\alpha}),0\}>0$, then they are removed from the set because they are infeasible. After the removal, a legitimate APF and approximate Pareto optimal set can be attained. The reader can refer to [69] for a similar technique that relates an SOP to an MOP.

V Numerical Examples

This section presents various system and control design examples using BMI approaches. Among the solution methods included for comparison, only BB methods involve global optimization. A detailed description of design problems and associated system parameters can be found in the appendices. Sections V-A and V-B examine feasibility problems and SOPs constrained by BMIs, respectively. Algorithm parameters $N_{nom}=40$, $N_{max}=160$, and $t_{max}=20$ were used, and 70 simulation runs were performed. For the BMI-based MOPs in Section V-C, the iteration number $t_{max}=300$ was used to produce APFs. These parameters were chosen based on a number of experiments in consideration of the algorithm convergence and computational time.

V-A Feasibility Problems with BMI Constraints

Table I presents our simulation results. See Appendix A for detailed problem descriptions. The “SR $\%$” represents the success rate of the proposed MRV solving these feasibility problems with BMI constraints. While the AM, ILMI, diffeomorphic state transformations, and two-step procedure were able to solve respective problems, our method successfully found solutions in a unified manner.

V-B SOPs with BMI Constraints

We compared the MRV with BB and path-following methods. Table II shows the numerical results in which the “Optimum, Mean,” and “Std” stand for the achieved optimal value, mean value, and standard deviation, respectively. In comparison with existing methods, the MRV yielded prominent improvement in problem SSS, and had similar results in problems LPVS and MCD.

To further assess the performance of the MRV, we used various models in COMPleib [6, 70], including aircraft models (AC), helicopter models (HE), reactor models (REA), decentralized interconnected systems (DIS), wind energy conversion models (WEC), terrain following models (TF), and academic test problems (NN). Spectral abscissa optimization problems, $H_{2}$ optimization problems, and $H_{\infty}$ optimization problems were investigated. The associated system under investigation has the following form:

[TABLE]

The closed-loop system of (32) using a static output feedback controller $\bm{u}=\bm{Fy}=\bm{FCx}$ can be written as

[TABLE]

The spectral abscissa optimization (or minimization) associated with (33) is formulated as [10]

[TABLE]

where

[TABLE]

is the spectral abscissa of $\bm{A}_{\bm{F}}$, $eig\{\bm{A}_{\bm{F}}\}$ represents the set of eigenvalues of $\bm{A}_{\bm{F}}$, $Re(\lambda)$ is the real part of $\lambda$, and matrix $\bm{F}$ represents the controller gain that must be determined. Because the objective function in (34) is neither smooth nor Lipschitz continuous, (34) is conventionally transformed into the BMI problem [10, 6]:

[TABLE]

where $\beta$ is related to the decay rate of the system.

To use our methodology, we compared (34) to (3), and let $\bm{\alpha}=\bm{F}$, $\bm{X}=\emptyset$ (no internal variable is involved), and $\bm{\mathcal{F}}(\bm{\alpha})=\alpha_{o}(\bm{A}_{\bm{F}})$. Table III presents the resulting performance of various solution methods.444In Tables III–V, the numerical results of existing methods HIFOO, LMIRank, PENBMI, and CCDM come from [4], and the results of ICAM are from [5]. The notation $\alpha_{o}(\bm{A})$ represents the spectral abscissa of $\bm{A}$. Minimization problems are considered. Values in the columns of Table III labeled with “Min” and “Mean” present the best possible performance and average performance during the simulation trials, respectively. The letter “x” means that no solution is found. Our approach performed excellently in approximately $73\%$ of test problems (marked in bold numbers) and yielded similar levels of performance in the remaining problems as compared with existing solution methods.

Remark 3

In Tables III–VI, the “Min” and “Mean” serve as performance metrics in different situations. Having the minimum “Min” in the results of the MRV implies that the proposed method outperforms existing solution methods in the best-case scenario. The best-case scenario can be related to the situation in which the computational complexity is not a concern. The best solution can then be obtained by a series of evaluations. This situation occurs in certain off-line applications and the value of “Min” can serve as a performance metric. By contrast, having the minimum “Mean” in the results of the MRV indicates that the proposed method is better than existing solution methods in average. When computational resources are limited, e.g., in certain online applications, the value of “Mean” can serve as a performance metric.

For $H_{2}$ and $H_{\infty}$ optimization, we use

[TABLE]

to represent the open- and closed-loop systems, respectively. The controller gain $\bm{F}$ was designed so that the $H_{2}$ norm of the closed-loop system, denoted by $||G_{c\ell}(\bm{F})||_{2}$, or the $H_{\infty}$ norm of the closed-loop system, denoted by $||G_{c\ell}(\bm{F})||_{\infty}$, was minimized while certain BMI constraints were satisfied (see Appendix B). Tables IV and V present the respective results.555In the tables, $||\bm{G}_{o\ell}||_{2}$ and $||\bm{G}_{o\ell}||_{\infty}$ represent the $H_{2}$ and $H_{\infty}$ norms of the open-loop system $\bm{G}_{o\ell}$, respectively. The notation “Inf” stands for “infinity.” The MRV outperformed existing solution methods in approximately $27.5\%$ and $47.8\%$ of test problems for $H_{2}$ and $H_{\infty}$ optimization, respectively (marked in bold numbers). For the remaining problems, it yielded similar levels of performance.

As testified in [4] and [5], the CCDM and ICAM were relatively robust compared with other existing solution methods, but they failed or made little progress towards a local solution in problems AC18, DIS5, PAS, and NN12 in COMPleib. By contrast, the MRV was able to find solutions in these problems, as shown in Table VI.

V-C MOPs with BMI Constraints

This subsection examines the ability of the MRV to produce APFs for controller designs involving multiple objectives. The first problem is the sparse linear constant output-feedback design described in (10) in which the matrices $\bm{A},\bm{B},$ and $\bm{C}$ are defined in [2] and [4]. The goal is to maximize the decay rate $\beta$ and minimize the entry values of the controller gain $\bm{F}$.

The second problem is an MO version of a mixed $H_{2}/H_{\infty}$ control problem (derived from (35) in [4]):

[TABLE]

where $\bm{A}_{\bm{F}}=\bm{A}+\bm{BFC}$, $\bm{C}_{\bm{F}}^{\bm{z}_{i}}=\bm{C}^{\bm{z}_{i}}+\bm{FC},i=1,2$, and the matrices $\bm{A},\bm{B},\bm{B}_{1},\bm{C}_{\bm{F}}^{\bm{z}_{1}}$ and $\bm{C}_{\bm{F}}^{\bm{z}_{2}}$ are defined in [2] and [4]. The $H_{2}$ and $H_{\infty}$ performance are related to $trace(\bm{Z})$ and $\gamma$, respectively. To apply the classification principle, we introduced a slack variable $\eta$ and imposed an additional constraint $trace(\bm{Z})\leq\eta^{2}$. The mixed $H_{2}/H_{\infty}$ design problem in (37) can then be transformed into

[TABLE]

According to the classification principle, $\bm{P}_{1},\bm{P}_{2}$, and $\bm{Z}$ must be included in the internal variable $\bm{X}$ because they have a constraint on positive definiteness ($\bm{Z}$ is located in a diagonal block of a positive-definite matrix and hence, it is positive-definite); since $\eta$ and $\gamma$ appear in the objective function, they are included in the external variable $\bm{\alpha}$; and finally, $\bm{F}$ must be included in $\bm{\alpha}$ as well so that for a fixed $\bm{\alpha}$, $\mathcal{BMI}(\bm{\alpha},\bm{X})<0$ becomes an LMI in the variable $\bm{X}$. Referring to the notations in (3), we let $\bm{\alpha}=(\eta,\gamma,\bm{F})$, $\bm{X}=(\bm{P}_{1},\bm{P}_{2},\bm{Z})$, and $\bm{\mathcal{F}}(\bm{\alpha})=[\eta\;\gamma]^{T}$.

Fig. 2 shows APFs obtained by solving (10) and (38). The APF in Fig. 2(a) is bent, implying that objectives are not heavily dependent. By choosing a design that corresponds to a vector in the knee region of the APF, it is possible to simultaneously improve both objectives, i.e., maximizing the decay rate and minimizing the values of entries of the controller gain. By contrast, the line shape of the APF in Fig. 2(b) indicates that simultaneous improvement in both objectives cannot be attained. Therefore, for the design problem in (38), we must sacrifice the $H_{\infty}$ performance to improve the $H_{2}$ performance, or vice versa. These two examples illustrate the ability of the MRV to produce APFs for MOPs constrained by BMIs. Since an APF can contain useful information about the relationships among objectives, applying the proposed methodology to solve BMI-constrained MOPs can be advantageous.

Remark 4

We illustrated how beneficial it can be by using the proposed methodology, consisting of variable classification, problem transformation, and algorithm integration, to solve BMI-constrained problems in system and control designs. The comparisons with existing BMI solution methods validated the effectiveness of the proposed methodology. A proof of concept was thus provided. It is worth noting that existing MOEAs cannot be directly applied to obtain solutions to our BMI-constrained problems such as (10) and (38) because they are not expressed in a standard form of an MOP. (The problem in (16) derived from the Equivalence Theorem is in the standard form.) Even if our variable classification and problem transformation have been performed so that the standard form has been obtained, there is little chance that existing MOEAs can solve the resulting MOPs because the external variable $\bm{\alpha}$ often contains the controller gain $\bm{F}$ that is still related to some matrix constraint. The proposed hybrid algorithm can have the power to produce solutions mainly because we have incorporated a pole-placement technique into the search engine (described in Section III-B) and designed a mechanism that ensures legitimate pole placement (described in Section IV-A).

VI Conclusion

In this paper, we proposed a solution method termed the MRV for BMI problems in system and control designs. By using this method, the associated decision variables are classified into external and internal variables according to the variable classification principle; the BMI problem is then transformed into an unconstrained optimization problem that has fewer decision variables; and finally, a hybrid algorithm termed HMOIA is applied to solve the unconstrained problem, yielding a feasible point, a solution, or a set of approximate Pareto optimal solutions depending on the dimension of the objective function. In our simulations, we compared the proposed MRV to various BMI solution methods and found that the MRV yielded excellent levels of performance in many benchmark problems, validating the proposed methodology. In contrast with some existing BMI solution methods, the MRV possesses the following advantages: it expresses decision variables in a vector form, which is convenient for controller designs; it avoids much effort such as problem reformulation or prior derivations, which can be heuristic and cumbersome; it performs global optimization instead of local search, which is essential because BMI problems are non-convex and have multiple local optima; and it can address multiple objectives simultaneously.

Appendix A Feasibility Problems with BMI Constraints

If not specified, the following equation numbers are those in the respective references.

ST (Stability Test, Sec. V-A of [11] with $\mu=0.1$): The problem was described in (5) of this paper, and the bounds $\tau_{\ell ij}\in[0,10]$ were used.

SIP (Stabilization of Inverted Pendulum, Sec. V-B of [11] with $\mu=0.001$ and BMIs in (13)):

[TABLE]

where $\tau_{\ell ij}\geq 0$, $\bm{F}_{i}$, and $\bm{P}_{i}>0$ are the decision variables. The external and internal variables were classified as $\bm{\alpha}=(\tau_{112},\tau_{121},\tau_{212},\tau_{221},\bm{F}_{1},\bm{F}_{2})$ and $\bm{X}=(\bm{P}_{1},\bm{P}_{2})$, respectively. The bounds $\tau_{\ell ij}\in[0,10]$ and $[\bm{F}_{i}]_{mn}\in[-10,10]$ were used.

SAFS-I (Stabilization of an Affine Fuzzy System, Sec. V of [7] with BMIs in (16.1) and (16.2)):

[TABLE]

and

[TABLE]

for $(i,j)=(1,1),(3,3),(1,2),(2,3)$, where

[TABLE]

The $\tau_{ij}\geq 0$, $\bm{F}_{i}$, and $\bm{P}>0$ are the decision variables. The external and internal variables were classified as $\bm{\alpha}=(\tau_{11},\tau_{33},\tau_{12},\tau_{23},\bm{F}_{1},\bm{F}_{2},\bm{F}_{3})$ and $\bm{X}=\bm{P}$, respectively. The bounds $\tau_{ij}\in[0,5]$ and $[\bm{F}_{i}]_{mn}\in[-5,5]$ were used.

SAFS-II (Stabilization of an Affine Fuzzy System, Sec. V of  [8] with BMIs in (14a) and (14b)):

[TABLE]

and

[TABLE]

for $(i,j)=(1,1),(3,1)$. The $\tau_{ij}\geq 0$, $\bm{F}_{i}$, and $\bm{P}>0$ are the decision variables. The external and internal variables were classified as $\bm{\alpha}=(\tau_{11},\tau_{31},\sigma_{1},\sigma_{3},\bm{F}_{1},\bm{F}_{2},\bm{F}_{3})$ and $\bm{X}=\bm{P}$, respectively. The bounds $\tau_{ij}\in[0,5],\sigma_{i}\in[-5,5]$ and $[\bm{F}_{i}]_{mn}\in[-5,5]$ were used.

OCS (Observer-based Control System, Sec. IV-A of [9] with BMIs in (16)):

[TABLE]

for $i=1,2,3,4$. The $\bm{F}_{i},\bm{G}_{i},i=1,2,3,4$, and $\bm{P}_{i}>0$, $i=1,2,$ are the decision variables, and the observer gains $\bm{L}_{i},i=1,2,3,4,$ are recovered by $\bm{L}_{i}=\bm{P}_{2}^{-1}\bm{G}_{i}$. The external and internal variables were classified as $\bm{\alpha}=(\bm{F}_{1},\bm{F}_{2},\bm{F}_{3},\bm{F}_{4})$ and $\bm{X}=(\bm{P}_{1},\bm{P}_{2},\bm{G}_{1},\bm{G}_{2},\bm{G}_{3},\bm{G}_{4})$, respectively. The bounds

[TABLE]

for $i=1,2,3,4$, were used.

Appendix B SOPs and MOPs with BMI Constraints

If not specified, the following equation numbers are those in the respective references.

LPVS (Linear Parameter-varying System, (28)–(32) in Sec. V-A of [29]): The problem was presented in (8) of this paper, and the bounds $\delta_{i}\in[0,1]$ and $\varsigma\in[0,10]$ were used.

SSS (Simultaneous State-feedback Stabilization, Sec. 4.2 [2]): A stabilizing state-feedback gain $\bm{F}$ exists if the optimum of

[TABLE]

is positive. We let $\bm{\alpha}=(\gamma_{1},\gamma_{2},\gamma_{3},\bm{F})$, $\bm{X}=(\bm{P}_{1},\bm{P}_{2},\bm{P}_{3})$, and $\bm{\mathcal{F}}(\bm{\alpha})=-\min\{\gamma_{1},\gamma_{2},\gamma_{3}\}$. The bounds $[\bm{F}]_{mn}\in[-50,50]$ and $\gamma_{i}\in[0,5]$ were used.

MCD (Mixed $H_{2}$/$H_{\infty}$ Controller Design, Sec. III-A of [15] with BMIs in (3)): We let $\bm{\alpha}=(\eta,\bm{K})$ and $\bm{X}=(\bm{P}_{1},\bm{P}_{2})$. The bounds $\eta\in[0,2]$ and $[\bm{K}]_{n}\in[-5,5]$ were used.

Spectral Abscissa Optimization: The bounds $[\bm{F}]_{mn}\in[-50,50]$ and $eig\{\bm{A}_{\bm{F}}\}\in\{\sigma+j\omega:(\sigma,\omega)\in[-20,0]\times[-20,20]\}$ were used, $S=3$ subspaces were adopted, and $\kappa_{s}$ was chosen from $\{1,0.5,0.1\}$ uniformly at random when each $\bm{\alpha}_{i}$ in (25) of this paper was constructed. The same setting was used in the $H_{2}$ optimization, $H_{\infty}$ optimization, and MOPs as well.

$H_{2}$ Optimization: We let $\bm{D}_{11}=0$ in $G_{c\ell}(\bm{F})$ and solved

[TABLE]

The term $\bm{B}_{1}\bm{B}_{1}^{T}$ was replaced by $\bm{B}_{1}\bm{B}_{1}^{T}+10^{-5}\bm{I}$ if it was not positive definite. We let $\bm{\alpha}=\bm{F}$, $\bm{X}=(\bm{Q},\bm{Y})$, and $\bm{\mathcal{F}}(\bm{\alpha})=||G_{c\ell}(\bm{F})||_{2}$ according to the classification principle. Given the value of the external variable $\bm{\alpha}=\bm{F}$, existing deterministic algorithms can be applied to evaluate $\bm{\mathcal{F}}(\bm{\alpha})=||G_{c\ell}(\bm{F})||_{2}$, e.g., the MATLAB routine $\mathrm{norm(syst,p)}$ with $\mathrm{syst}=G_{c\ell}(\bm{F})$ and $\mathrm{p}=2$ can be used. To facilitate numerical comparisons and expedite the solving process, we used the following settings for both $H_{2}$ and $H_{\infty}$ optimization: the objective values $\tilde{\bm{\mathcal{F}}}:=[10^{5}\;10^{5}+\alpha_{o}(\bm{A}_{\bm{F}})]^{T}$ were assigned without further evaluation of $\lambda^{*}(\bm{\alpha})$ whenever $\alpha_{o}(\bm{A}_{\bm{F}})\geq 0$; and the bounds $eig\{\bm{A}_{\bm{F}}\}\in\{\sigma+j\omega:(\sigma,\omega)\in[-20,0]\times[-20,20]\}$ were not used (only bounds on $[\bm{F}]_{mn}$ were used) when problem NN11 in COMPleib was solved.

$H_{\infty}$ Optimization: We solved

[TABLE]

We let $\bm{\alpha}=\bm{F}$, $\bm{X}=(\gamma,\bm{Y})$, and $\bm{\mathcal{F}}(\bm{\alpha})=||G_{c\ell}(\bm{F})||_{\infty}$. Given the value of the external variable $\bm{\alpha}=\bm{F}$, the value of $\bm{\mathcal{F}}(\bm{\alpha})=||G_{c\ell}(\bm{F})||_{\infty}$ can be determined using deterministic algorithms, e.g., the MATLAB routine $\mathrm{norm(syst,p)}$ with $\mathrm{syst}=G_{c\ell}(\bm{F})$ and $\mathrm{p=inf}$ can be used.

MOPs: We used the bound $\beta\in[0,1.5]$ in the sparse linear constant output-feedback design, and the bounds $\eta\in[0,2],\gamma\in[1,5]$ in the mixed $H_{2}$/$H_{\infty}$ design.