Multimodeling Control via System Balancing

Abstract

A new approach in multimodeling strategy is proposed. Multimodel strategies in which control agents use different simplified models of the same system are being developed using balancing transformation and the corresponding order reduction concepts. Traditionally, the multimodeling concept was studied using the ideas of multitime scales (singular perturbations) and weak subsystem coupling. For all reduced-order models obtained, a Linear Quadratic Gaussian (LQG) control problem was solved. Different order reduction techniques were compared based on the values of the optimized criteria for the closed-loop case where the full-order balanced model utilizes regulators calculated to be the optimal for various reduced-order models. The results obtained were demonstrated on a real-world example: a multiarea power system consisting of two identical areas, that is, two identical power plants.

1. Introduction

Large-scale systems have been the subject of research work for several decades [1–20]. Some order reduction techniques were developed for the singularly perturbed class of systems, based on different mathematical procedures, such as graph metric [21]. The concept of multimodel strategies for large-scale systems is originated from [17]. According to that concept, a large-scale system may be controlled by several independent agents using several simplified models of the system. In the singular perturbation [17], methodology has been used to develop multimodel strategies by exploiting the nature of the system that has two fast subsystems mutually weakly coupled and both strongly connected to the slow subsystem. The basic contribution of [17] has been to establish a set of conditions under which a multimodel strategy is well posed in the sense that the performance of the control system is close to the performance that would have been obtained had the control strategy been designed by a single controller knowing the exact model of the overall system. Multimodel strategies [17] have been studied and employed in several papers, either within the context of singular perturbations or in different mathematical set-ups [14]. Conditions for a multimodel strategy to be well-posed are investigated in [14] for a linear quadratic Gaussian (LQG) optimal control problem [22].

The multimodeling structure used in the classical approach via singular perturbations and weak coupling is defined by a linear dynamic system that has one slow and two fast subsystems [11]. The fast subsystems are strongly connected to the slow subsystem and weakly connected (or not connected) among themselves. This structure describes well the dynamics of several real-world systems, for example, power systems [15, 17] and automobiles [2]. The corresponding multimodeling representation [11, 17] is defined by (1.1): 𝑥̇₀𝜀(𝑡)₁𝑥̇₁𝜀(𝑡)₂𝑥̇₂=𝐴(𝑡)⁰⁰𝐴₀₁𝐴₀₂𝐴₁₀𝐴₁₁𝜀₃𝐴₁₂𝐴₂₀𝜀₃𝐴₂₁𝐴₂₂𝑥₀𝑥(𝑡)₁𝑥(𝑡)₂+𝐵(𝑡)⁰¹𝐵₀₂𝐵₁₁𝜀₃𝐵₁₂𝜀₃𝐵₂₁𝐵₂₂𝑢₁𝑢(𝑡)₂𝑥(𝑡),𝑦(𝑡)=𝐶𝑥⁰(𝑡)₁𝑥(𝑡)₂𝑢(𝑡)+𝐷¹𝑢(𝑡)₂,(𝑡) (1.1)(1.2) where 𝑥₀∈ℜ^𝑛0 are slow state variables, 𝑥₁∈ℜ^𝑛1, 𝑥₂∈ℜ^𝑛2 are fast state variables, 𝑢₁∈ℜ^𝑚1, 𝑢₂∈ℜ^𝑚2 are control inputs, 𝑦∈ℜ^𝑝 are outputs, 𝜀₁, 𝜀₂ are small positive singular perturbation parameters, and 𝜀₃ is small weak-coupling parameter. (For 𝜀₃=0, (1.1) describes multiparameter singularly perturbed system (MSPS) studied in detail in literature [5–11].)

For the purpose of deterministic optimal control of the above multimodeling structure, the quadratic performance criterion has to be minimized by the proper choice of the control variables 𝑢₁(𝑡) and 𝑢₂(t). The performance criterion for the linear quadratic Gaussian optimization [22] is given by 1𝐽=2₀^+∞[𝑥^𝑇(𝑡)𝑄𝑥(𝑡)+𝑢^𝑇(𝑡)𝑅𝑢(𝑡)]𝑑𝑡,𝑄=𝑄^𝑇≥0,𝑅=𝑅^𝑇>0, (1.3) where 𝑥𝑥(𝑡)=𝑥⁰(𝑡)₁𝑥(𝑡)₂𝑢(𝑡),𝑢(𝑡)=¹𝑢(𝑡)₂𝑄(𝑡),𝑄=⁰⁰𝑄₀₁𝑄₀₂𝑄^𝑇₀₁𝑄₁₁0𝑄^𝑇₀₂0𝑄₂₂𝑅,𝑅=¹00𝑅₂,𝑄=𝑞^𝑇𝑞𝑞=⁰¹𝑞₁₁0𝑞₀₂0𝑞₂₂^𝑇𝑞₀₁𝑞₁₁0𝑞₀₂0𝑞₂₂=𝑞^𝑇₀₁𝑞₀₁+𝑞^𝑇₀₂𝑞₀₂𝑞^𝑇₀₁𝑞₁₁𝑞^𝑇₀₂𝑞₂₂𝑞^𝑇₁₁𝑞₀₁𝑞^𝑇₁₁𝑞₁₁0𝑞^𝑇₂₂𝑞₀₂0𝑞^𝑇₂₂𝑞₂₂. (1.4) In the general multimodeling case, all zero-elements in matrices R and Q can be replaced by O(ε) elements. (𝑂(𝜀^𝑖) is defined by 𝑂(𝜀^𝑖)<𝑐𝜀^𝑖, where 𝑐 is a bounded constant, 𝑖 is a real number, and 𝜀=‖[𝜀₁𝜀₂]‖.)

To solve the multimodeling problem one proceeds with constructing two different models of (1.1), obtained by setting 𝜀₁=0, which leads to the first model for the first controller, and by setting 𝜀₂=0, which produces the second model for the second controller. In order to simplify equations, without loss of generality, small coupling parameter 𝜀₃ is set to zero.

The fast dynamics of the other subsystem is approximated by an algebraic equation (the corresponding 𝜀_𝑖 is set to zero).

Two approximations of the original model derived from (1.1) are 𝑥̇₀𝜀(𝑡)₁𝑥̇₁0=𝐴(𝑡)⁰⁰𝐴₀₁𝐴₀₂𝐴₁₀𝐴₁₁0𝐴₂₀0𝐴₂₂𝑥₀𝑥(𝑡)₁𝑥(𝑡)₂+𝐵(𝑡)⁰¹𝐵₀₂𝐵₁₁00𝐵₂₂𝑢₁𝑢(𝑡)₂,𝑥̇(𝑡)₀0𝜀(𝑡)₂𝑥̇₂=𝐴(𝑡)⁰⁰𝐴₀₁𝐴₀₂𝐴₁₀𝐴₁₁0𝐴₂₀0𝐴₂₂𝑥₀𝑥(𝑡)₁𝑥(𝑡)₂+𝐵(𝑡)⁰¹𝐵₀₂𝐵₁₁00𝐵₂₂𝑢₁𝑢(𝑡)₂.(𝑡) (1.5) The above equations can be rewritten as: 𝑥̇₀𝜀(𝑡)₁𝑥̇₁=(𝑡)(𝐴₀₀−𝐴₀₂𝐴⁻¹₂₂𝐴₂₀)𝐴₀₁𝐴₁₀𝐴₁₁𝑥₀𝑥(𝑡)₁+𝐵(𝑡)⁰¹(𝐵₀₂−𝐴₀₂𝐴⁻¹₂₂𝐵₂₂)𝐵₁₁0𝑢₁𝑢(𝑡)₂,𝑥̇(𝑡)₀𝜀(𝑡)₂𝑥̇₂=(𝑡)(𝐴₀₀−𝐴₀₁𝐴⁻¹₁₁𝐴₁₀)𝐴₀₂𝐴₂₀𝐴₂₂𝑥₀𝑥(𝑡)₂+(𝑡)(𝐵₀₁−𝐴₀₁𝐴⁻¹₁₁𝐵₁₁)𝐵₀₂0𝐵₂₂𝑢₁𝑢(𝑡)₂,(𝑡) (1.6) leading to two different models of the original system. The algebraic equations defined in (1.5) are used in (1.3) by each controller to form their own performance criterion. Such simplified criteria are optimized by each controller via the corresponding reduced-order model and the obtained control strategies form the multimodeling strategy. The multimodeling strategy is well posed if it is 𝑂(𝜀) close to the global optimal control strategy obtained by performing direct optimization on the original system and the original performance criterion, as shown in [14].

2. The Use of Balancing Transformation for the System Order Reduction

Robust order reduction based on the use of balancing transformation has been described in [20]. Concisely, consideration is given to a linear, time invariant system as in (2.1): 𝑑𝑥(𝑡)𝑑𝑡=𝐴𝑥(𝑡)+𝐵𝑢(𝑡),𝑥(𝑡₀)=𝑥₀,𝑦(𝑡)=𝐶𝑥(𝑡)+𝐷𝑢(𝑡), (2.1) where 𝑥(𝑡) is an n-dimensional state vector, 𝑢(𝑡) is an m-dimensional input vector, and 𝑦(𝑡) is a p-dimensional output vector.

For a linear, time invariant system (2.1) the corresponding transfer function for the open-loop system is given by (2.2): 𝐺(𝑠)=𝐶(𝑠𝐼−𝐴)⁻¹𝐵+𝐷. (2.2)

It is assumed that the system (2.1) is asymptotically stable and that a G(s) is of minimal realization.

Assumption. A system is asymptotically stable, a pair (𝐴,𝐵) is controllable, and a pair (𝐴,𝐶) is observable.

The controllability and observability Gramians of the original system (2.1) satisfy the algebraic equations of Lyapunov [15, 16]: 𝑃𝐴^𝑇+𝐴𝑃+𝐵𝐵^𝑇=0,𝑄𝐴+𝐴^𝑇𝑄+𝐶^𝑇𝐶=0. (2.3)

For a system that is controllable and observable, both controllability Gramian and observability Gramian are positive definite matrices, 𝑃>0, 𝑄>0.

The balancing transformation is applied on the space vector in order to achieve that the controllability and the observability Gramians are identical and diagonal, that is, 𝑥_𝑏(𝑡)=𝑇𝑥(𝑡),det(𝑇)≠0⇒𝑑𝑥_𝑏(𝑡)𝑑𝑡=𝐴_𝑏𝑥_𝑏(𝑡)+𝐵_𝑏𝑦𝑢(𝑡),_𝑏(𝑡)=𝐶_𝑏𝑥_𝑏(𝑡)+𝐷_𝑏𝐴𝑢(𝑡)=𝑦(𝑡),_𝑏=𝑇𝐴𝑇⁻¹,𝐵_𝑏=𝑇𝐵,𝐶_𝑏=𝐶𝑇⁻¹,𝐷_𝑏𝑃=𝐷,_𝑏=𝑄_𝑏=Σ=diag{𝜎₁,𝜎₂,…,𝜎_𝑛},𝜎₁≥𝜎₂≥⋯≥𝜎_𝑛>0, (2.4)(2.5)(2.6) where 𝜎_𝑖 are known as the Hankel singular values (HSVs).

Assuming that the original system is controllable and observable, a balanced system will also be both controllable and observable, since the similarity transformation preserves controllability and observability of the system that was transformed [15, 16]. Hence all 𝜎_𝑖 are positive. Furthermore, both original and balanced systems are of minimal realization. The transfer function of the balanced system given by 𝐺_𝑏(𝑠)=𝐶_𝑏(𝑠𝐼−𝐴_𝑏)⁻¹𝐵_𝑏+𝐷=𝐺(𝑠) (2.7) stays unchanged thanks to a coordinate change through a nonsingular transformation. The balanced controllability and observability Gramians satisfy the following algebraic Lyapunov equations: Σ𝐴^𝑇_𝑏+𝐴_𝑏Σ+𝐵_𝑏𝐵^𝑇_𝑏=0⇔𝑃_𝑏𝐴^𝑇_𝑏+𝐴_𝑏𝑃_𝑏+𝐵_𝑏𝐵^𝑇_𝑏=0,Σ𝐴_𝑏+𝐴^𝑇_𝑏Σ+𝐶^𝑇_𝑏𝐶_𝑏=0⇔𝑄_𝑏𝐴_𝑏+𝐴^𝑇_𝑏𝑄_𝑏+𝐶^𝑇_𝑏𝐶_𝑏=0. (2.8)

The idea of the order reduction through balancing transformation can be linked with the canonical system decomposition. It was shown that the system’s modes that were either uncontrollable or unobservable did not appear in the system transfer function. In [15, 16] it is shown that the system modes that are both weakly controllable and weakly observable have little influence on the system dynamics; so they can be neglected.

However, it was noticed that those modes which are weakly controllable and well observable cannot be neglected as can neither be neglected those modes that are well controllable and weakly observable. Let us assume that the balanced system (2.4)–(2.6) is partitioned in the following way: 𝐴_𝑏=𝐴₁₁𝐴₁₂𝐴₂₁𝐴₂₂,𝐵_𝑏=𝐵₁₁𝐵₂₂,𝐶_𝑏=𝐶₁₁𝐶₂₂,𝐷_𝑏Σ=𝐷,Σ=¹00Σ₂,Σ₁=diag{𝜎₁,𝜎₂,…,𝜎_𝑟},Σ₂=diag{𝜎_𝑟+1,𝜎_𝑟+2,…,𝜎_𝑛}. (2.9) Assuming that 𝜎_𝑟>𝜎_𝑟+1, balanced truncation produces a system of lower order, r, defined by 𝑑𝑥₁(𝑡)𝑑𝑡=𝐴₁₁𝑥₁(𝑡)+𝐵₁₁𝑢(𝑡),𝑦(𝑡)=𝐶₁₁𝑥₁(𝑡)+𝐷𝑢(𝑡), (2.10) and the corresponding transfer function of the reduced-order system is 𝐺₁₁(𝑠)=𝐶₁₁(𝑠𝐼−𝐴₁₁)⁻¹𝐵₁₁+𝐷. (2.11)

The reduced-order system attained in this way is both controllable and observable since all corresponding HSVs are positive. Furthermore, the reduced-order system is balanced and asymptotically stable. It was shown in literature, for example, [20] that the 𝐻_∞ norm for the reduced-order system, obtained through the truncation procedure given above, satisfies the condition: ‖𝐺(𝑠)−𝐺₁₁(𝑠)‖_∞≤2(𝜎_𝑟+1+𝜎_𝑟+2+⋯+𝜎_𝑛). (2.12)

It was noticed that the reduced-order system obtained through the balanced truncation procedure gives very good approximation of the original system in the case of pulse input for both control signals (good spectra approximation on higher frequencies) but shows a considerable, steady-state error in the case of step input (poor spectra approximation on lower frequencies). This error is due to the fact that the original system and the reduced order system have different DC gains. Actually, after the above-described truncation through balancing transformation, most of the spectra on lower frequencies are kept and also some of the spectra on higher frequencies, but some of the spectra on lower frequencies are lost as well as most of the spectra on higher frequencies. By eliminating part of the spectra on lower frequencies (which occur in the neglected part of the system—state variables 𝑥₂(𝑡)) we have produced gain that differs from the gain of the original system that was balanced. This discrepancy was eliminated in [16] where a technique of balanced residualization was proposed that produced an accurate (exact) DC gain and a very good spectra approximation on lower frequencies and sometimes even on middle frequencies. It should be noted that in [20] a residualisation technique was also used. An improved truncation method that preserves the exact DC gain value as in the original system is given in [16] and applied in [19].

3. Multimodeling via System Balancing

Two multimodeling order reduction concepts are mentioned above: the first one uses the small-parameter idea and the second one is based on balancing transformation, presented through two specific methods—balanced truncation and balanced residualization. Here, we will present the idea of multimodeling via system balancing, which is more general than multimodeling via singular perturbation since the latter requires special structures of the original model. Multimodeling via singular perturbation is performed on the assumption that two fast weakly connected subsystems, whose system matrices are invertible, are both strongly coupled to the slow subsystem while no such assumption is necessary for the multimodeling via system balancing.

Let us consider a time invariant linear system represented by (3.1): 𝑥̇₀𝑥̇(𝑡)₁𝑥̇(𝑡)₂𝑥(𝑡)=𝐴𝑥⁰(𝑡)₁𝑥(𝑡)₂𝑢(𝑡)+𝐵¹𝑢(𝑡)₂,𝑥(𝑡)𝑦(𝑡)=𝐶𝑥⁰(𝑡)₁𝑥(𝑡)₂𝑢(𝑡)+𝐷¹𝑢(𝑡)₂,(𝑡) (3.1) where 𝐴𝐴=⁰⁰𝐴₀₁𝐴₀₂𝐴₁₀𝐴₁₁𝐴₁₂𝐴₂₀𝐴₂₁𝐴₂₂𝐵,𝐵=⁰¹𝐵₀₂𝐵₁₁𝐵₁₂𝐵₂₁𝐵₂₂. (3.2) The balanced model is derived as in (2.4)-(2.5), that is, 𝑥_𝑏(𝑡)=𝑇𝑥(𝑡),det(𝑇)≠0⇒𝑑𝑥_𝑏(𝑡)𝑑𝑡=𝐴_𝑏𝑥_𝑏(𝑡)+𝐵_𝑏𝑦𝑢(𝑡),_𝑏(𝑡)=𝐶_𝑏𝑥_𝑏(𝑡)+𝐷_𝑏𝐴𝑢(𝑡)=𝑦(𝑡),_𝑏=𝑇𝐴𝑇⁻¹,𝐵_𝑏=𝑇𝐵,𝐶_𝑏=𝐶𝑇⁻¹,𝐷_𝑏𝑃=𝐷,_𝑏=𝑄_𝑏=Σ=diag{𝜎₁,𝜎₂,…,𝜎_𝑛},𝜎₁≥𝜎₂≥⋯≥𝜎_𝑛>0. (3.3) The HSVs are the same as the HSVs of the original model (1.1) and are sorted in a descending order so that they can help in determining how to partition the model: the first few state space variables for the slow subsystem and the following two for the two fast parts, as in (3.4).

The balanced model will be partitioned in the following manner (3.4)-(3.5): 𝑥̇_𝑏0𝑥̇(𝑡)_𝑏1𝑥̇(𝑡)_𝑏2=𝐴(𝑡)^𝑏00𝐴_𝑏01𝐴_𝑏02𝐴_𝑏10𝐴_𝑏11𝐴_𝑏12𝐴_𝑏20𝐴_𝑏21𝐴_𝑏22𝑥_𝑏0𝑥(𝑡)_𝑏1𝑥(𝑡)_𝑏2+𝐵(𝑡)^𝑏01𝐵_𝑏02𝐵_𝑏11𝐵_𝑏12𝐵_𝑏21𝐵_𝑏22𝑢₁𝑢(𝑡)₂𝑦(𝑡)_𝑏,(𝑡)=𝐶_𝑏0𝐶_𝑏1𝐶_𝑏2𝑥_𝑏0𝑥(𝑡)_𝑏1𝑥(𝑡)_𝑏2+𝐷(𝑡)^𝑏1𝐷_𝑏2𝑢₁𝑢(𝑡)₂(𝑡)=𝑦(𝑡), (3.4)(3.5) where 𝑥_𝑏0∈ℜ^𝑛0 are slow or the common-core state variables of the balanced model; 𝑥_𝑏1∈ℜ^𝑛1, 𝑥_𝑏2∈ℜ^𝑛2 are state variables of the balanced model, having a small amount of total energy (𝑛₀+𝑛₁+𝑛₂=𝑛); 𝑢₁∈ℜ^𝑚1, 𝑢₂∈ℜ^𝑚2 are control inputs of the balanced and original model (𝑚₁+𝑚₂=𝑚); 𝑦_𝑏∈ℜ^𝑝, 𝑦∈ℜ^𝑝 are outputs of the balanced and original model; respectively, and 𝐴_𝑏𝑖𝑗,𝐵_𝑏𝑖𝑘,𝐶_𝑏𝑗, and 𝐷_𝑏𝑘 (𝑖,𝑗=0,1,2 and 𝑘=1,2) are submatrices of the balanced system matrices 𝐴_𝑏,𝐵_𝑏,𝐶_𝑏, and 𝐷_𝑏, having the corresponding dimensions.

Instead of zero-submatrices in [12] a more general case is considered here, where submatrices 𝐴_𝑏12, 𝐴_𝑏21, 𝐵_𝑏12, and 𝐵_𝑏21 are nonzero matrices, as it is usually the case in the real-world systems.

The corresponding quadratic performance criterion which has to be minimized is 𝐽_𝑏=12₀^+∞[𝑥^𝑇_𝑏(𝑡)𝑄_𝑏𝑥_𝑏(𝑡)+𝑢^𝑇(𝑡)𝑅_𝑏𝑢(𝑡)]𝑑𝑡,𝑄_𝑏=𝑄^𝑇_𝑏≥0,𝑅_𝑏=𝑅^𝑇_𝑏>0, (3.6) where, without loss of generality, we assume that 𝑄_𝑏 and 𝑅_𝑏 are block-diagonal matrices: 𝑄_𝑏=𝑄_𝑏00000𝑄_𝑏11000𝑄_𝑏22,𝑅_𝑏=𝑅_𝑏1100𝑅_𝑏22. (3.7)

Now the decomposition of the balanced model (3.4) is performed into two reduced-order models, both of them composed of the common core—a subsystem corresponding to slow modes of the balanced model, and one of the two different remaining subsystems, corresponding to the modes having a small amount of total energy, namely, the first reduced-order model and the second reduced-order model. Cross-coupling matrices 𝐴_𝑏12,𝐴_𝑏21 will be neglected though they have nonzero values. The first reduced-order model obtained is 𝑥̇_𝑏0𝑥̇(𝑡)_𝑏1(𝑡)=𝐴_1𝑟𝑥_𝑏0𝑥(𝑡)_𝑏1(𝑡)+𝐵_1𝑟𝑢₁𝑢(𝑡)₂,𝑦(𝑡)_𝑏1(𝑡)=𝐶_𝑏0𝐶_𝑏1𝑥_𝑏0𝑥(𝑡)_𝑏1(𝑡)+𝐷_𝑏1𝑢₁𝐴(𝑡),_1𝑟=𝐴_𝑏00𝐴_𝑏01𝐴_𝑏10𝐴_𝑏11,𝐵_1𝑟=𝐵_𝑏01𝐵_𝑏02𝐵_𝑏11𝐵_𝑏12. (3.8)(3.9)The second reduced-order model is 𝑥̇_𝑏0𝑥̇(𝑡)_𝑏2(𝑡)=𝐴_2𝑟𝑥_𝑏0𝑥(𝑡)_𝑏2(𝑡)+𝐵_2𝑟𝑢₁𝑢(𝑡)₂,𝑦(𝑡)_𝑏2(𝑡)=𝐶_𝑏0𝐶_𝑏2𝑥_𝑏0𝑥(𝑡)_𝑏2(𝑡)+𝐷_𝑏2𝑢₂(𝑡), (3.10) where 𝐴_2𝑟=𝐴_𝑏00𝐴_𝑏02𝐴_𝑏20𝐴_𝑏22,𝐵_2𝑟=𝐵_𝑏01𝐵_𝑏02𝐵_𝑏12𝐵_𝑏22. (3.11)

It is interesting to note that the first and the second reduced-order models in (3.8) and (3.10) both have HSVs corresponding to the common core modes of the original system and two complement sets of the HSVs corresponding to the modes of the original system that have a small amount of total energy, if we keep all parts of the 𝐵_𝑏 matrix.

If, however, coupling submatrices 𝐵_𝑏12 and 𝐵_𝑏21 are set to zero-matrices of the corresponding dimensions, as suggested in [11], the HSVs obtained for the first and the second subsystems will differ from those of the original full-order model, and the reduced-order models will be as in (3.12)–(3.13) and (3.14). These two types of reduced-order models would still contain all control signals of the original and balanced models.

The first reduced-order model is obtained as in [11, 17]:𝑥̇_𝑏0𝑥̇(𝑡)_𝑏1(𝑡)=𝐴_1𝑟𝑥_𝑏0𝑥(𝑡)_𝑏1(𝑡)+𝐵_1𝑟𝑢₁𝑢(𝑡)₂,𝑦(𝑡)_𝑏1(𝑡)=𝐶_𝑏0𝐶_𝑏1𝑥_𝑏0𝑥(𝑡)_𝑏1(𝑡)+𝐷_𝑏1𝑢₁(𝑡), (3.12) where 𝐴_1𝑟=𝐴_𝑏00𝐴_𝑏01𝐴_𝑏10𝐴_𝑏11,𝐵_1𝑟=𝐵_𝑏01𝐵_𝑏02𝐵_𝑏110. (3.13)

The second reduced-order model is obtained as in [11, 17]:𝑥̇_𝑏0𝑥̇(𝑡)_𝑏2(𝑡)=𝐴_2𝑟𝑥_𝑏0𝑥(𝑡)_𝑏2(𝑡)+𝐵_2𝑟𝑢₁𝑢(𝑡)₂,𝑦(𝑡)_𝑏2(𝑡)=𝐶_𝑏0𝐶_𝑏2𝑥_𝑏0𝑥(𝑡)_𝑏2(𝑡)+𝐷_𝑏2𝑢₂𝐴(𝑡)._2𝑟=𝐴_𝑏00𝐴_𝑏02𝐴_𝑏20𝐴_𝑏22,𝐵_2𝑟=𝐵_𝑏01𝐵_𝑏020𝑏_𝑏22. (3.14)

The same approximation is done for the performance criterion (1.3); hence two performance criteria are needed, which leads to multicriteria optimization problem. Depending on the actual problem setup, very often described by differential games, the two controllers find their own optimal strategies and apply such strategies to the global system defined by (1.1). In such a way, the multimodeling strategy is well posed if the performance criterion under the multimodeling strategy is 𝑂(𝜎_𝑟+1/𝜎_𝑟) close to the global optimal control strategy obtained by performing direct optimization of the original performance criterion for the original system.

The following two criteria to be optimized now are 𝐽₁=12₀^+∞{𝑥_𝑏0(𝑡)𝑥_𝑏1(𝑡)𝑄₁𝑥_𝑏0𝑥(𝑡)_𝑏1+𝑢(𝑡)¹(𝑡)𝑢₂(𝑡)𝑅₁𝑢₁𝑢(𝑡)₂𝑄(𝑡)}𝑑𝑡,₁=𝑄^𝑇₁≥0,𝑅₁=𝑅^𝑇₁𝐽>0,₂=12₀^+∞{𝑥_𝑏0(𝑡)𝑥_𝑏2(𝑡)𝑄₂𝑥_𝑏0𝑥(𝑡)_𝑏2+𝑢(𝑡)¹(𝑡)𝑢₂(𝑡)𝑅₂𝑢₁𝑢(𝑡)₂𝑄(𝑡)}𝑑𝑡,₂=𝑄^𝑇₂≥0,𝑅₂=𝑅^𝑇₂>0. (3.15)(3.16)𝑄₁ and 𝑄₂ are block-diagonal matrices and 𝑅₁ and 𝑅₂ are submatrices of 𝑅_𝑏 having corresponding dimensions: 𝑄₁=𝑄^𝑇₁=𝑄_𝑏0000𝑄_𝑏11≥0,𝑅₁=𝑅^𝑇₁=𝑅_𝑏11𝑄>0,₂=𝑄^𝑇₂=𝑄_𝑏0000𝑄_𝑏22≥0,𝑅₂=𝑅^𝑇₂=𝑅_𝑏22>0. (3.17)

The optimal control problem we refer to, following the suggestions in [2, 14, 17], is to minimize a convex sum of 𝐽₁ and 𝐽₂ given in (3.15) and (3.16), that is, for some 𝛾₁ and 𝛾₂𝐽_𝑀𝑀𝐶=𝛾₁𝐽₁+𝛾₂𝐽₂,𝛾₁+𝛾₂=1,𝛾₁>0,𝛾₂>0. (3.18) This corresponds to a Pareto optimal cooperative strategy [6].

In Pareto optimal strategy a situation is considered in which decision makers should decide on their strategies through mutual cooperation [14, 17]. The essence of this is that no variation from a Pareto optimal strategy can decrease the costs of either of the decision makers. Let each decision maker have a quadratic cost functional as in (3.15) and (3.16). A Pareto solution is a pair 𝑢_opt1(𝑡),𝑢_opt2(𝑡) which minimizes (3.18) for some 𝛾₁ and 𝛾₂ [6]. The optimal feedback solution to (3.15) and (3.16) is given by (3.19): 𝑢_opt𝑢(𝑡)=^opt1𝑢(𝑡)_opt2(𝑡)=−𝑅_𝑏⁻¹𝐵^𝑇_𝑏𝑃_𝑏𝑥_𝑏(𝑡). (3.19) Here 𝑃_𝑏 is the positive semidefinite stabilizing solution of the algebraic Riccati equation: 𝐴^𝑇_𝑏𝑃_𝑏+𝑃_𝑏𝐴_𝑏+𝑄_𝑏−𝑃_𝑏𝑆_𝑏𝑃_𝑏=0,𝑆_𝑏=𝐵_𝑏𝑅_𝑏⁻¹𝐵^𝑇_𝑏. (3.20)

The optimal state regulator is a special case of this problem where the decision makers agree on a choice of 𝛾₁ and 𝛾₂ as weighting factors. Without prejudice to the generality, we chose 𝛾₁=𝛾₂=0.5. We expect that the approximation of the optimization criteria would be 12(𝐽₁+𝐽₂)≈𝐽_𝑏, (3.21) and according to [17] it is to be expected that: 𝐽_𝑀𝑀𝐶=12𝐽₁+12𝐽₂𝜎+𝑂(^𝑟+1𝜎_𝑟)≈𝐽. (3.22)

The original result is obtained in [9], where instead of 𝑂(𝜎_𝑟+1/𝜎_𝑟), 𝑂(‖𝜇‖²) was used, where 𝜇 stood for a norm of [𝜀₁𝜀₂𝜀₃], in specific ‖𝜇‖=𝜀²₁+𝜀²₂+𝜀²₃. Still adhering to this idea but including criteria for balancing truncation, instead of 𝑂(‖𝜇‖²) we can use 𝑂(𝜎_𝑟+1/𝜎_𝑟) having in mind that the order r is chosen where 𝜎_𝑟+1 is of order of percent of the 𝜎_𝑟. In [6] a recursive algorithm is developed for solving multiparameter Riccati equations with the rate of convergence 𝑂(‖𝜇‖^𝑖) and the rate of accuracy for the near-optimal strategy is 𝑂(‖𝜇‖^𝑖+1) where i is the iteration number.

The required solution of the algebraic Riccati equation (3.20) is based on the standard assumption [5, 12].

Assumption. The triples (𝐴_1𝑟,𝐵_1𝑟,chol(𝑄₁)) and (𝐴_2𝑟,𝐵_2𝑟,chol(𝑄₂)) are stabilizable-detectable. (Here chol(Q) is the Cholesky decomposition of a matrix 𝑄.)

Assumption. The pairs (𝐴_1𝑟,chol(𝑄₁)) and (𝐴_2𝑟,chol(𝑄₂)) are detectable-observable.

The matrices (𝐴_𝑖𝑟,𝐵_𝑖𝑟), 𝑖=1,2 are given by (3.23): 𝐴_1𝑟=𝐴_𝑏00𝐴_𝑏01𝐴_𝑏10𝐴_𝑏11,𝐵_1𝑟=𝐵_𝑏01𝐵_𝑏02𝐵_𝑏110,𝐴_2𝑟=𝐴_𝑏00𝐴_𝑏02𝐴_𝑏20𝐴_𝑏22,𝐵_2𝑟=𝐵_𝑏01𝐵_𝑏020𝐵_𝑏22. (3.23) The optimal feedback solution to (3.12)–(3.13) and (3.15) is given by 𝑢_thefirstopt𝑢(𝑡)=^{firstopt1(𝑡)}𝑢_{firstopt2(𝑡)}=−𝑅₁⁻¹𝐵^𝑇_1𝑟𝑃₁𝑥_𝑏0𝑥(𝑡)_𝑏1(𝑡), (3.24) where 𝑃₁ is the positive semidefinite stabilizing solution of the corresponding algebraic Riccati equation (3.25): 𝐴^𝑇_1𝑟𝑃₁+𝑃₁𝐴_1𝑟+𝑄₁−𝑃₁𝑆₁𝑃₁=0,𝑆₁=𝐵_1𝑟𝑅₁⁻¹𝐵^𝑇_1𝑟. (3.25) The optimal feedback solution to (3.14) and (3.16) is given by 𝑢_thesecondopt𝑢(𝑡)=^secondopt1𝑢(𝑡)_secondopt2(𝑡)=−𝑅₂⁻¹𝐵^𝑇_2𝑟𝑃₂𝑥_𝑏0𝑥(𝑡)_𝑏2(𝑡), (3.26) where 𝑃₂ is the positive semidefinite stabilizing solution of the corresponding algebraic Riccati equation: 𝐴^𝑇_2𝑟𝑃₂+𝑃₂𝐴_2𝑟+𝑄₂−𝑃₂𝑆₂𝑃₂=0,𝑆₂=𝐵_2𝑟𝑅₂⁻¹𝐵^𝑇_2𝑟. (3.27)

It is intuitively clear that a further reduction could be made for some models. If the control signals in the original model are weakly coupled, then it is possible to form two subsystems as described above, however the first one having only inputs 𝑢₁ and the second subsystem having only inputs 𝑢₂. This would mean that the neglected modes having a small amount of total energy are considered to have reached their steady-state values, rather than as changing variables. Models of the type (1.1) are sound candidates for this type of reduction. This would be as if the submatrices 𝐵_𝑏02 in (3.13) and 𝐵_𝑏01 in (3.14) were set to zero-matrices of the corresponding dimensions, resulting in (3.28)–(3.31) and (3.32)–(3.35), respectively. The first reduced-order model obtained by neglecting the control variable 𝑢₂ in (3.12)–(3.13) is 𝑥̇_𝑏0𝑥̇(𝑡)_𝑏1=𝐴(𝑡)^𝑏00𝐴_𝑏01𝐴_𝑏10𝐴_𝑏11𝑥_𝑏0𝑥(𝑡)_𝑏1+𝐵(𝑡)^𝑏010𝐵_𝑏110𝑢₁𝑢(𝑡)₂(𝑡), (3.28) or 𝑥̇_𝑏0𝑥̇(𝑡)_𝑏1(𝑡)=𝐴_1𝑟𝑥_𝑏0𝑥(𝑡)_𝑏1(𝑡)+𝐵_1𝑟𝑢₁𝑦(𝑡),_𝑏1(𝑡)=𝐶_𝑏0𝐶_𝑏1𝑥_𝑏0𝑥(𝑡)_𝑏1(𝑡)+𝐷_𝑏1𝑢₁𝐴(𝑡),_1𝑟=𝐴_𝑏00𝐴_𝑏01𝐴_𝑏10𝐴_𝑏11,𝐵_1𝑟=𝐵_𝑏01𝐵_𝑏11. (3.29)(3.30)(3.31)

The second reduced-order model obtained by neglecting the control variables 𝑢₁ in (3.14) is presented in (3.32)–(3.35): 𝑥̇_𝑏0𝑥̇(𝑡)_𝑏2=𝐴(𝑡)^𝑏00𝐴_𝑏02𝐴_𝑏20𝐴_𝑏22𝑥_𝑏0𝑥(𝑡)_𝑏2+(𝑡)0𝐵_𝑏020𝐵_𝑏22𝑢₁𝑢(𝑡)₂(𝑡), (3.32) or 𝑥̇_𝑏0𝑥̇(𝑡)_𝑏2(𝑡)=𝐴_2𝑟𝑥_𝑏0𝑥(𝑡)_𝑏2(𝑡)+𝐵_2𝑟𝑢₂𝑦(𝑡),_𝑏2(𝑡)=𝐶_𝑏0𝐶_𝑏2𝑥_𝑏0𝑥(𝑡)_𝑏2(𝑡)+𝐷_𝑏2𝑢₂𝐴(𝑡),_2𝑟=𝐴_𝑏00𝐴_𝑏02𝐴_𝑏20𝐴_𝑏22,𝐵_2𝑟=𝐵_𝑏02𝐵_𝑏22. (3.33)(3.34)(3.35)

The same approximation is done for the performance criterion (3.15)-(3.16); hence two performance criteria are obtained, as in (3.36)-(3.37) and (3.38)-(3.39), each of them having its own control signal. Taking into consideration multicriteria optimization, both of these criteria would be regarded as having equal importance. Again, the multimodeling strategy will be considered to be well posed if the performance criterion under the multimodeling strategy is O(𝜎_𝑟+1/𝜎_𝑟) close to the global optimal control strategy obtained by performing direct optimization on the original system and the original performance criterion.

The two criteria to be optimized now are 𝐽_thefirst=12₀^+∞{𝑥_𝑏0(𝑡)𝑥_𝑏1(𝑡)𝑄₁𝑥_𝑏0𝑥(𝑡)_𝑏1(𝑡)+𝑢₁(𝑡)^𝑇𝑅₁𝑢₁𝑄(𝑡)}𝑑𝑡,₁=𝑄^𝑇₁=𝑄_𝑏0000𝑄_𝑏11≥0,𝑅₁=𝑅^𝑇₁=𝑅_𝑏11𝐽>0,_thesecond=12₀^+∞{𝑥_𝑏0(𝑡)𝑥_𝑏2(𝑡)𝑄₂𝑥_𝑏0𝑥(𝑡)_𝑏2(𝑡)+𝑢₂(𝑡)^𝑇𝑅₂𝑢₂𝑄(𝑡)}𝑑𝑡,₂=𝑄^𝑇₂=𝑄_𝑏0000𝑄_𝑏22≥0,𝑅₂=𝑅^𝑇₂=𝑅_𝑏22>0. (3.36)(3.37)(3.38)(3.39)

In the simulation example it has been chosen that 𝑄₁ and 𝑄₂ are unity matrices of the corresponding dimensions, while 𝑅₁ and 𝑅₂ are unity matrices of the corresponding dimensions multiplied by some weighing factor which is in accordance with the suggestions made in [17]: 𝑄₁=𝑄^𝑇₁=𝐼_𝑛0+𝑛1≥0,𝑅₁=𝑅^𝑇₁=Const⋅𝐼_𝑚1𝑄>0,₂=𝑄^𝑇₂=𝐼_𝑛0+𝑛2≥0,𝑅₂=𝑅^𝑇₂=Const⋅𝐼_𝑚2>0. (3.40)(3.41) The optimal feedback solution to (3.36)-(3.37) and (3.40) is given by (3.42): 𝑢_thefirstopt(𝑡)=−𝑅₁⁻¹𝐵^𝑇_1𝑟𝑃_thefirst𝑥_𝑏0𝑥(𝑡)_𝑏1(𝑡),𝐵_1𝑟=𝐵_𝑏01𝐵_𝑏11, (3.42) where 𝑃_thefirst is the positive semidefinite stabilizing solution of the algebraic Riccati equation: 𝐴^𝑇_1𝑟𝑃_thefirst+𝑃_thefirst𝐴_1𝑟+𝑄₁−𝑃_thefirst𝑆_thefirst𝑃_thefirst=0,𝑆_thefirst=𝐵_1𝑟𝑅₁⁻¹𝐵^𝑇_1𝑟, (3.43) while the optimal feedback solution to (3.38)-(3.39) and (3.41) is given by: 𝑢_thesecondopt(𝑡)=−𝑅₂⁻¹𝐵^𝑇_2𝑟𝑃_thesecond𝑥_𝑏0𝑥(𝑡)_𝑏2(𝑡),𝐵_2𝑟=𝐵_𝑏02𝐵_𝑏22. (3.44) Here 𝑃_thesecond is the positive semidefinite stabilizing solution of the algebraic Riccati equation: 𝐴^𝑇_2𝑟𝑃_thesecond+𝑃_thesecond𝐴_2𝑟+𝑄₂−𝑃_thesecond𝑆_thesecond𝑃_thesecond=0,𝑆_thesecond=𝐵_2𝑟𝑅₂⁻¹𝐵^𝑇_2𝑟. (3.45) It could be explored at a later stage what the sufficient and necessary conditions should be for the existence of solutions to Riccati equations (3.43) and (3.45). Some useful ideas and efficient iterative solutions with the existence conditions could be found in the work of Mukaidani et al. [5–11].

The next step is to implement optimal control from (3.42) and (3.44) to the full-order model. Optimal regulators (3.42) and (3.44) are easier to design than an optimal regulator for the full-order model (3.6). Since both control agents are of equal importance to the full-order model, approximation of the optimal control strategy could be 𝑢_approxopt𝑢(𝑡)=^thefirstopt𝑢(𝑡)_thesecondopt=(𝑡)−𝑅₁⁻¹𝐵^𝑇_1𝑟𝑃_thefirst𝑥_𝑏0𝑥(𝑡)_𝑏1(𝑡)−𝑅₂⁻¹𝐵^𝑇_2𝑟𝑃_thesecond𝑥_𝑏0𝑥(𝑡)_𝑏2(𝑡), (3.46) as was shown in the simulation example.

The optimal control problem we refer to, following the suggestions in [2, 14, 17], is to minimize a convex sum of 𝐽_thefirst and 𝐽_thesecond given in (3.36) and (3.38) for some 𝛾₁ and 𝛾₂: 𝐽_approximated=𝛾₁𝐽_thefirst+𝛾₂𝐽_thesecond,𝛾₁+𝛾₂=1,𝛾₁>0,𝛾₂>0. (3.47) This corresponds to a Pareto optimal cooperative strategy [6].

In Pareto optimal strategy we have considered a situation in which decision makers decide on their strategies through mutual cooperation [14, 17]. The essence of this approach is that no variation from a Pareto optimal strategy can decrease the costs of any of the decision makers. Let each decision maker have a quadratic cost functional as in (3.36) and (3.38). A Pareto solution is a pair 𝑢_thefirstopt(𝑡),𝑢_thesecondopt(𝑡) which minimized (3.47) for some 𝛾₁ and 𝛾₂ [6]. The optimal state regulator is a special aspect of this problem where the decision makers agree on the choice of 𝛾₁ and 𝛾₂ as weighting factors.

Without prejudice to the generality of the above considerations, we chose 𝛾₁=𝛾₂=0.5. According to [17] it is to be expected that the approximation of the optimization criteria would be 𝐽_approximated=𝐽_bal.IandII=12(𝐽_thefirst+𝐽_thesecond)≈𝐽_𝑏. (3.48)

4. Example

The methods for the order reduction displayed here were tested on an example known from literature and taken from [17] with some modification. Such example is a state space model of a power system consisting of two interconnected identical areas, where each area consists of one plant. The model behavior was simulated in the open-loop as well as in the closed-loop.

For the model example a system was chosen having two inputs and one output. Based on HSVs a decision was made on what should be taken for the reduced order while the order of the slow or common-core subsystem and of the two fast parts (i.e., parts having a small amount of total energy) was chosen as suggested in [17]. Several reduced-order models were produced, using the methods as mentioned and described above.

The efficiency of these approximations was compared first in the open-loop case for typical input functions: impulse, step, ramp, and sine. In the open-loop case a comparison was also made with the balanced model for all frequency characteristics (both magnitude and phase spectra) of all available transfer functions. Control strategy in this paper was different from the one used in [17]. Optimal linear quadratic regulators (LQRs) were designed for original and balanced models as well as for the ones of reduced order. Gain matrices of the reduced-order models were expanded by zeros to obtain the full order, and these reduced-order regulators were used to close the loop on the balanced model. Different approximations were compared with respect to the values of optimized LQG criteria.

The model in the state-space representation is of the order 𝑛=9 with the astatism of the second order therefore being unstable and is described by means of the matrices:,𝐵𝐴=004.50100000004.5−1000000−0.050−0.10.1000000−0.050.1000.100032.7−32.70000000000−550000−40000−10000000000−55000−400000−10_𝑇=.000.100010000000.1000010 (4.1)

Inspired by [1, page 227, Case 3] and [17] we have chosen 𝐶=110001010,𝐷=00. (4.2)

For the model chosen for simulation, state-space variables as in [17] are as follows: 𝑥₁ and 𝑣₁ the integral of the area control error for the area 1, 𝑥₂ and 𝑣₂ the integral of the area control error for the area 2, 𝑥₃ and Δ𝑓₁ frequency variation for the area 1, 𝑥₄ and Δ𝑓₂ frequency variation for the area 2, 𝑥₅ and Δ𝑃₁₂ tie-line power flow variation, 𝑥₆ and Δ𝑃_𝐺1 turbine output variation for the plant 1 or the area 1, 𝑥₇ and Δ𝑎₁ turbine valve position variation for the plant 1 or the area 1, 𝑥₈ and Δ𝑃_𝐺2 turbine output variation for the plant 2 or the area 2, and 𝑥₉ and Δ𝑎₂ turbine valve position variation for the plant 2 or the area 2, and control signals are as follows: 𝑢₁ and Δ𝑃_𝑐1 speed changer variation for the plant 1 or the area 1, and 𝑢₂ and Δ𝑃_𝑐1 speed changer variation for the plant 2 or the area 2.

It is significant to note that the astatism is not inherent to the system—it is induced by including the integrals of the area control error (ACE) 𝑣₁ and 𝑣₂ into state vector. So the simulation could be performed with the modification of this model where these two variables would be omitted from the state space vector but retained in the output or outputs, resulting in the model having full-order 𝑛=7.

Here we have chosen another modification. To make the model in the simulation example stabilized, matrix 𝐴 is changed according to the prescribed degree of stability: 𝐴=𝐴+𝛼𝐼₉, (4.3) where 𝛼 is chosen to be −10 and 𝐼₉ denotes unity matrix of dimension 9. Other system matrices remained unchanged. Matrix 𝐴 which is changed to have the prescribed degree of stability is 𝐴=−1004.50100000−1004.5−1000000−10.050−0.10.1000000−10.050.1000.100032.7−32.7−10000000000−1550000−40000−20000000000−155000−400000−20. (4.4)

However, the system obtained in this way is still unobservable—only four of its modes are observable. From the matrix dimensions in (4.4) it is obvious that 𝑛=9 is the order of the system.

Table 1 shows HSVs of the original and balanced models. On the basis of the HSVs a decision was made that the reduced order could be (and so it was chosen to be) 𝑟=2, since the third Hankel singular value is more than 400 (410.4444) times smaller than the second one. For the reduced-order model of order 𝑟=2 good quality approximation will be achieved, both in the open- and in closed-loop cases. However, in this manner we would omit state space variables that are important to be retained.

Table 1 

Hankel singular values.

Multimodeling allows us to reduce the order of the model while keeping all state space variables accounted for, only decoupled into two intersecting subsets.

As for the multimodeling, the reduced-order system which consisted of a common core subsystem and subsystem corresponding to modes having a small amount of total energy was called the first while the second referred to the one which consisted of a common core subsystem and subsystem corresponding to the remaining modes, having a smaller amount of total energy.

According to the suggestions made in [17] the slow or common core subsystem order was chosen to be of order 𝑛₀=5, and the orders of the two fast subsystems, that is, subsystems corresponding to modes with a small amount of total energy were 𝑛₁=2, and 𝑛₂=2, respectively; so the first and the second subsystems were of the same order, 𝑛₀+𝑛₁=𝑛₀+𝑛₂=7. For the first and the second reduced-order models a multimodeling was performed through balanced truncation. In this way two subsystems were obtained having one input, with u₁ being input to the first subsystem, as in (3.28)–(3.31), and 𝑢₂ being input to the second one, as in (3.32)–(3.35).

Table 1 contains also the HSVs for the first and the second reduced-order models (ROMs) obtained through multimodeling using balanced truncation. It seems that the largest four HSVs for the first subsystem and the second one are identical. However, the difference between them is of the order 10⁻¹⁵ or less.

Each of the model approximations as well as the original one was tested on the open-loop case for typical input functions: impulse, step, ramp, and sine. In the time domain all of the approximations produce good and similar performances, except for both the first and the second subsystems obtained through multimodeling, that exhibit a considerable degree of difference in behavior when compared to models in the cases of ramp and sine input signals. However, the behavior of the first model and the second one is almost identical.

The optimal gain matrices were subsequently calculated for the original model, its balanced equivalent, and the four reduced-order models: one obtained through balanced truncation, one obtained through balanced residualization, whereas two reduced-order models were comprised of the slow subsystem and one fast subsystem—called the first and the second—consisted of the slow subsystem and the rest of the fast subsystem. Matrices needed for optimization criteria, chosen in accordance with those from [17], and used for the original and balanced nonreduced model and for various reduced-order models are in (4.5): 𝑄=𝐼₉,𝑅=20𝐼₂,𝑄_truncated=𝑄_rezidualized=𝐼₂,𝑅_truncated=𝑅_rezidualized=20𝐼₂,𝑄_thefirst=𝑄_thesecond=𝐼₇,𝑅_thefirst=𝑅_thesecond=20. (4.5) Here 𝐼_𝑘 denotes unity matrix of dimension 𝑘.

Simulations were carried out using MATLAB. Proportional regulator gain matrices were computed for the original and the balanced models of the full-order and the truncated and residualized reduced-order models of order 2 as well as for the first and the second reduced models of order 7.

Corresponding matrices for the first and the second subsystems, 𝐾_thefirst and 𝐾_thesecond, are close to submatrices of 𝐾_balanced and of the necessary dimensions:𝐾_balanced=−3.28⋅10⁻³−1.561⋅10⁻³0.34⋅10⁻⁴−0.15⋅10⁻⁴−0.33⋅10⁻³0.89⋅10⁻³−0.32⋅10⁻⁵−0.111⋅10⁻³−2.65⋅10⁻³−3.28⋅10⁻³−1.561⋅10⁻³0.34⋅10⁻⁴−0.15⋅10⁻⁴0.586⋅10⁻³−0.89⋅10⁻³0.32⋅10⁻⁵0.111⋅10⁻³2.65⋅10⁻³,𝐾_truncated=−3.2778⋅10⁻³−1.5597⋅10⁻³−3.2778⋅10⁻³−1.5597⋅10⁻³𝐾_rezidualized=−3.2785⋅10⁻³−1.5635⋅10⁻³−3.2785⋅10⁻³−1.5635⋅10⁻³,𝐾_thefirst=−3.28⋅10⁻³−1.561⋅10⁻³0.34⋅10⁻⁴−0.15⋅10⁻⁴−0.397⋅10⁻³−0.78⋅10⁻⁴−0.247⋅10⁻³,𝐾_thesecond=−3.280⋅10⁻³−1.561⋅10⁻³0.34⋅10⁻⁴−0.15⋅10⁻⁴0.586⋅10⁻³−0.111⋅10⁻³2.234⋅10⁻³. (4.6)

Finally, a regulator is composed of the 𝐾_thefirst and 𝐾_thesecond as in 𝐾_bal.IandII=−3.28⋅10⁻³−1.561⋅10⁻³0.34⋅10⁻⁴−0.15⋅10⁻⁴−0.397⋅10⁻³−0.78⋅10⁻⁴−0.247⋅10⁻³00−3.28⋅10⁻³−1.561⋅10⁻³0.34⋅10⁻⁴−0.15⋅10⁻⁴0.586⋅10⁻³00−0.111⋅10⁻³2.234⋅10⁻³. (4.7)

The closed-loop systems were also tested for various input functions: impulse, step, ramp and sine. Regulators, calculated on reduced order models were subsequently used to form the closed loop on the balanced full-order model.

The corresponding Lyapunov equations were solved and the traces of their solutions were compared to derive which of the closed-loop models obtained in this manner is the closest to the optimal case with a regulator having gain matrix 𝐾_balanced. The values of the optimal criteria for the balanced model with its optimal 𝑃 regulator and the values of suboptimal criteria for the balanced model with an optimal regulator for the various reduced-order models, are seemingly identical, as in Criteria𝐵: Criteria𝐵=0.45840.45840.45840.45840.45840.4584. (4.8) From the values of Criteria𝐵 it appears that we have reached the desired approximation of the optimization criteria: 𝐽_approximated=𝐽_bal.IandII=12𝐽_thefirst+12𝐽_thesecond≈0.4584. (4.9) Absolute errors with respect to balanced full-order optimal case areCriteriaError𝐵=00.768⋅10⁻⁵0.772⋅10⁻⁵0.1402⋅10⁻⁴0.1030⋅10⁻⁴0.469⋅10⁻⁵. (4.10)

However, the difference in the optimal criteria value is in the CriteriaError𝐵 and for the balanced model with optimal regulator for the truncated model it is 0.768⋅10^-5.

For the balanced model with optimal regulator for the residualized the difference in the optimal criteria value is 0.772⋅10^-5. The difference in the optimal criteria value for the balanced model with optimal regulator for the first subsystem is 0.1402⋅10^-4. For the balanced model with optimal regulator for the second subsystem the difference is 0.1030⋅10^-4.

For the regulator composed of the 𝐾_thefirst and 𝐾_thesecond gain matrices as described above, the difference between the corresponding suboptimal criteria and the optimal one is the smallest of all tested reduced-order models: 0.469⋅10^-5. This concurs with the quality approximation expected from the results of [17].

The corresponding relative errors with respect to the optimal criteria value for the balanced system are CriteriaRelativeError𝐵=0%0.0017%0.0017%0.0031%0.0022%0.0010%. (4.11) Apparently, a combination of the first and the second optimal gain matrices produces performance of the closed-loop model that is the closest to the balanced full-order model. The relative difference in the optimal criteria value is in the Criteria Relative Error B and, respectively, for the balanced model using regulator optimal for the truncated model it is 0.0017%. As for the balanced model with optimal regulator for the residualized, it is identical to the above value, that is, 0.0017%. For the balanced model using regulator optimal for the first subsystem relative criteria error is 0.0031%, and finally, for the balanced model with optimal regulator for the second subsystem it is 0.0022%.

For the regulator composed of the 𝐾_thefirst and 𝐾_thesecond as described above, the relative difference between the corresponding suboptimal criteria and the optimal one is the smallest of all tested reduced-order models: 0.0010%. This result confirms that the approximation achieved in this way is of a good quality.

5. Conclusion

Many order reduction techniques have been developed throughout the past few decades, and this problem will stay under consideration for as long as the engineering practice keeps developing. Here are performed and combined two already established techniques for order reduction by means of multimodeling and balancing. The novel approach has given very good results in the open-loop as well as in the closed-loop model approximation for real engineering purposes. In those cases where physical interpretation of the state space variables requires that not a single variable should be lost, multimodeling provides a way to reduce order and simplify the regulator design. The cases of this kind are the real-world models of power systems and cars, for example. The choice of order reduction technique can be made more appropriate provided that all kinds of information about the original system and the related constraints are incorporated.