*Article* **Symbolic Regulator Sets for a Weakly Nonlinear Discrete Control System with a Small Step**

**Yulia Danik \* and Mikhail Dmitriev**

Federal Research Center "Computer Science and Control" of the Russian Academy of Sciences, 119333 Mocow, Russia; mdmitriev@mail.ru

**\*** Correspondence: yuliadanik@gmail.com

**Abstract:** For a class of discrete weakly nonlinear state-dependent coefficient (SDC) control systems, a suboptimal synthesis is constructed over a finite interval with a large number of steps. A onepoint matrix Padé approximation (*PA*) of the solution of the initial problem for the discrete matrix Riccati equation is constructed based on the state-dependent Riccati equation (SDRE) approach and the asymptotics by the small-step of the boundary layer functions method. The symmetric gain coefficients matrix for Padé control synthesis is constructed based on the one-point *PA*. As a result, the parametric closed-loop control is obtained. The results of numerical experiments illustrate, in particular, the improved extrapolation properties of the constructed regulator, which makes the algorithm applicable in control systems for a wider range of parameter variation.

**Keywords:** discrete control systems; weakly nonlinear systems; small step; the SDRE approach; matrix discrete Riccati equation; the boundary layer functions method; Padé approximation; finite time interval

#### **Citation:** Danik, Y.; Dmitriev, M. Symbolic Regulator Sets for a Weakly Nonlinear Discrete Control System with a Small Step. *Mathematics* **2022**, *10*, 487. https://doi.org/10.3390/ math10030487

Academic Editors: Ioannis G. Stratis and Quanxin Zhu

Received: 9 November 2021 Accepted: 30 January 2022 Published: 2 February 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

#### **1. Introduction**

In the literature, much attention is paid to the construction of optimal control laws for nonlinear systems and the corresponding approximate methods for their calculation. This can be explained by several factors; on the one hand, the greater accuracy of the description of dynamic systems in applications leads to the increase of their mathematical models' dimension, and on the other hand, the calculations often need to be carried out in real time. This is especially true for finding feedback laws in nonlinear control systems, where the consideration of even weak nonlinearity in constructing synthesis laws based on linear control laws can lead to a significant improvement in the value of the quality criterion. For linear control systems, the Kalman algorithm is often used, which allows stabilizing feedbacks to be built that, in addition to the asymptotic stability of closed systems, provide the optimality by the quadratic quality criterion.

The application field of the Kalman algorithm has been expanded to nonlinear control systems by the so-called state-dependent Riccati equation (SDRE) approach for continuous (see [1–3]) and discrete (see [4–9]) cases, where the systems are formally represented as linear systems in terms of state and control, the coefficients of the matrices are the functions of the state vector (state-dependent coefficients (SDC) systems), and the quality criterion is quadratic, but the quadratic forms matrices in the criterion can also be state-dependent. In [1], the design procedure for the SDRE approach is described and its capabilities are illustrated on a benchmark problem. Reference [2] contains a detailed review on SDRE: the theory developed to date, characteristics, advantages, and open issues. In [4–6], the discrete state-dependent Riccati equation (D-SDRE) approach is described, and the corresponding control algorithms are presented and tested. In [3,7–9], the development of the SDRE approach is proposed based on the asymptotic theory.

Often, the mathematical models in applied control problems contain parameters, the variation of which generates a family of admissible controls and the corresponding trajectories, which leads to new problems of their approximate analytical study and description. An approximation is a procedure for choosing the optimal approximating function from a certain class of functions, that describes the behavior of the given function. Approximation allows the numerical characteristics and qualitative properties of an object to be studied, reducing the problem into studying simpler or more convenient objects. Moreover, the solution of these problems significantly reduces the time for rational control synthesis and selection. Various methods of approximation theory including splines and fractionalrational Padé approximations [10,11], where the latter are constructed using asymptotic expansions by the corresponding small parameters, can be used for such problems. In particular, in the works [10,11], the possibility of the construction of symbolic families of stabilizing regulators for some classes of regularly and singularly perturbed nonlinear continuous control systems are demonstrated by constructing Padé approximations (*PA*) of the gain coefficients matrices [10,11], associated with asymptotic expansions of the matrix algebraic Riccati equations solutions into regular series in integer powers of the corresponding small parameter, both for the case of one-point and two-point *PA*, where the latter uses not only the asymptotic expansions by integer powers of the parameter but also the asymptotic expansions over the inverse powers of a positive small parameter. A regulator is a closed-loop system, used to maintain the desired system output, usually, zero.

In this paper, discrete nonlinear control systems with a small step are considered within the framework of the SDRE approach for continuous systems on a finite time interval [12]. It is shown that an approximate symbolic description of the family of gain matrices in the feedback control loop using *PA* can be obtained if the step in the discrete system is sufficiently small. Note that the optimal control problem for a discrete-time system with a small step is singularly perturbed and that is why the zones of exponentially decreasing boundary layers can arise near the boundary points in the trajectories and in the solution of problems for finding gain coefficients in the feedback circuit. This happens because the degenerate solution, for zero value of the parameter, does not account for the boundary conditions and there is a large discrepancy from the exact solution near the boundary points. The boundary functions (members of the boundary series) are significant in the neighborhood of the boundary points and together with the members of the regular series, describe the behavior of the solution in the boundary layer, and outside the boundary layer, they rapidly decrease. The solution of the perturbed problem outside these boundary layer zones is in a small neighborhood of the solution of degenerate (limiting) problems, one of which, in this case, is the matrix algebraic Riccati equation of the corresponding discrete linear-quadratic optimal control problem [13], where the coefficients of all the matrices may be the functions of the state vector.

For the first time, the discrete problems with a small step were considered as singularly perturbed by A.B. Vasil'eva with the help of the boundary layer functions method (BLFM) (see [14,15]). The asymptotic methods and singular perturbations theory can also be successfully applied to control problems, for example, see reviews [16,17]. The discrete optimal control problem was first considered as singularly perturbed in [18] on the example of a linear-quadratic control problem with a small step, where its asymptotics by BLFM is constructed. The asymptotics of the solution of singularly perturbed discrete nonlinear optimal control problems with a small step was constructed by BLFM in [19], using the so-called direct scheme [20].

In this paper, the asymptotic solution of the corresponding initial singularly perturbed problem for the discrete matrix Riccati equation with coefficients weakly dependent on the state and the corresponding one-point *PA* regulator is constructed using the SDRE approach. Only one asymptotic expansion of the matrix Riccati equation solution by a small parameter, which equals the time step size and the multiplier in front of the non-linearities, is used to construct the *PA* for the gain matrix in the feedback loop. *PAs* are constructed based on asymptotic expansions and are actively used in applications since they often extend the range of the parameter variation, where *PAs* approach the exact solutions with given accuracy and reproduce the qualitative picture of the solutions. The use of a qualitative

picture allows the approximate solutions with higher accuracy to be obtained by using the results of the asymptotic analysis as an initial approximation for traditional nonlinear programming algorithms to minimize the residuals of the equations.

The results of numerical experiments are presented, which demonstrate the possibility of using this approach for the construction of stabilizing regulators for nonlinear discrete systems, and also show that the proposed algorithm is applicable for discrete control systems for a wider interval of the step value and the perturbation parameter variation.

The outline of the paper is as follows. In Section 2, we will describe the discrete state-dependent Riccati equation approach (SDRE) for small step discrete pseudo-linear control systems and construct the uniform second-order asymptotic approximation of the solution of the singularly perturbed initial problem for the difference Riccati equation using the boundary layer functions method (BLFM). In Section 3, we formulate the algorithm for the discrete one-point Padé regulator design. In Section 4, we will review the numerical experiments that demonstrate the work of the proposed control algorithms.

We start by listing definitions and notations: *<sup>R</sup>* := (−∞, <sup>∞</sup>), *<sup>R</sup>n*×*n*—vector space of n-by-n matrices, *<sup>T</sup>*—the transpose operation, ⊗—the Kronecker product, *En*—the identity matrix of size *n* × *n*, *P*—the regular series of *P*, Π*P*—the boundary series of *P*, *Pk*—*k*th term of the *P* series, the*L*<sup>2</sup> matrix norm.

#### **2. An SDRE Approach for Small Step Discrete Control Systems**

*2.1. Asymptotic Expansion of the Discrete Riccati Equation Solution*

Let's consider an affine (linear in control) discrete system,

$$\mathbf{x}(t+\varepsilon) = \overset{\frown}{A}(\mathbf{x}) + B(\mathbf{x})\boldsymbol{\mu}, \; \mathbf{x}(0) = \mathbf{x}\_{0\prime} \tag{1}$$

where *x*(*t*) is an n-dimensional state, *u*(*t*)—is an *r*-dimensional control, *<sup>t</sup>* <sup>∈</sup> *<sup>T</sup>* <sup>=</sup> {*<sup>t</sup>* : *<sup>t</sup>* <sup>=</sup> *<sup>k</sup>ε*, *<sup>k</sup>* <sup>=</sup> 0, 1, . . . , *<sup>N</sup>* <sup>−</sup> <sup>1</sup>} <sup>⊂</sup> {*<sup>t</sup>* : 0 <sup>≤</sup> *<sup>t</sup>* <sup>≤</sup> <sup>1</sup>}, *<sup>N</sup>* <sup>=</sup> *<sup>T</sup> <sup>ε</sup>* , *ε* > 0, is a small time step, used as a small parameter. Further, to transform the system into a formally linear form, we use the method of "extended linearization" (see [1,21,22]), which gives a non-unique representation in the vector case, by presenting *<sup>A</sup>*(*x*) in the form *A*(*x*) = *A*(*x*)*x*. By analogy, let us call the well-known heuristic technique for the introduction of a small parameter in the system right-hand side the "extended perturbation" method. In the last case, a small parameter is introduced in matrices *A*(*x*) and *B*(*x*) instead of a selected coefficient that approximately equals 1.

Let's demonstrate the last approach by transforming an arbitrary nonlinear function Θ(*x*), for example, as follows: <sup>Θ</sup>(*x*) = <sup>Θ</sup><sup>0</sup> <sup>+</sup> <sup>1</sup> · (Θ(*x*) <sup>−</sup> <sup>Θ</sup>0) = <sup>Θ</sup><sup>0</sup> <sup>+</sup> *<sup>ε</sup> <sup>δ</sup>* · (Θ(*x*) − Θ0) = Θ<sup>0</sup> + *ε*Θ1(*x*), where Θ1(*x*) = <sup>1</sup> *<sup>δ</sup>* ·(Θ(*x*) <sup>−</sup> <sup>Θ</sup>0) and 1 <sup>≡</sup> *<sup>ε</sup> <sup>δ</sup>* , but, at the same time, *δ* = *ε*, and now *δ* becomes a constant independent of *ε*. So, this technique by analogy with the previous one will be called the "extended" perturbation technique.

After application of these two approaches—"extended" linearization and "extended" perturbation—system (1) may be presented in the form

$$\begin{array}{c} \mathbf{x}(t+\varepsilon) = A(\mathbf{x}, \varepsilon)\mathbf{x} + B(\mathbf{x}, \varepsilon)u, \; \mathbf{x}(0) = \mathbf{x}\_{0\prime} \\ \mathbf{x}(t) \in X \in \mathbb{R}^n, \; u \in \mathbb{R}^r, \; t = k\varepsilon, k = 0, 1, \ldots, N - 1, \; 0 < \varepsilon \le \varepsilon\_{0\prime} \end{array} \tag{2}$$

where *A*(*x*,*ε*) = *A*<sup>0</sup> + *εA*1(*x*), *B*(*x*,*ε*) = *B*<sup>0</sup> + *εB*1(*x*), *ε*<sup>0</sup> is a small enough positive number, *<sup>A</sup>*0, *<sup>B</sup>*<sup>0</sup> are constant matrices, *<sup>A</sup>*0, *<sup>A</sup>*1(*x*) ∈ *<sup>R</sup>n*×*n*, *<sup>B</sup>*0, *<sup>B</sup>*1(*x*) ∈ *<sup>R</sup>n*×*<sup>r</sup>* , ∀*<sup>x</sup>* ∈ *<sup>X</sup>* ⊂ *<sup>R</sup><sup>n</sup>* and *<sup>X</sup>* is a certain fixed bounded state space subset. For sufficiently small values of *ε*, systems (2) are called singularly perturbed and pseudo-linear in state and control with regularly perturbed coefficients so that near the interval endpoints the solutions of (2) may have the boundary layers behavior.

Let's consider a cost function that measures the system performance for different controls in order to compare them and select the most favorable one;

$$J(\boldsymbol{\mu}) = \frac{1}{2} \mathbf{x}^T(\boldsymbol{N}) \mathbf{F} \mathbf{x}(\boldsymbol{N}) + \frac{1}{2} \sum\_{k=0}^{N-1} \left( \mathbf{x}^T(k\varepsilon) Q(\mathbf{x}, \varepsilon) \mathbf{x}(k\varepsilon) + \boldsymbol{\mu}^T(k\varepsilon) \mathbf{R} \boldsymbol{u}(k\varepsilon) \right) \to \min\_{\boldsymbol{\mu}}, \tag{3}$$

where *Q*(*x*,*ε*) = *Q*<sup>0</sup> + *εQ*1(*x*) + *ε*2*Q*2(*x*) > 0, *R* > 0, *Q*0, *R* are constant matrices, *Q*<sup>0</sup> > 0, *Q*1(*x*) > 0, *Q*2(*x*) > 0 ∀*x* ∈ *X*, *ε* ∈ (0,*ε*0) and *F* > 0.

To construct a feedback control for discrete systems on a finite time interval, we will use the necessary optimality conditions here, in contrast to [12], where a dynamic programming scheme was used in the continuous case.

For problems (1) and (2), we introduce the Hamiltonian

$$H(\mathbf{x}, u, \boldsymbol{\psi}, t) = \boldsymbol{\psi}^T(t)[A(\mathbf{x}, \varepsilon)\mathbf{x}(t) + B(\mathbf{x}, \varepsilon)u(t)] - \frac{1}{2}[\mathbf{x}^T(t)Q(\mathbf{x}, \varepsilon)\mathbf{x}(t) + u^T(t)R\_0u(t)]. \tag{4}$$

Using the necessary optimality conditions, we have

$$\mathbf{x}(t+\varepsilon) = A(\mathbf{x}, \varepsilon)\mathbf{x} + B(\mathbf{x}, \varepsilon)\mathbf{u},\tag{5}$$

$$\begin{split} \boldsymbol{\Psi}(t) &= \frac{\partial \boldsymbol{H}(\boldsymbol{x}(t), \boldsymbol{\upmu}(t), \boldsymbol{\upmu}(t+\varepsilon), t)}{\partial \boldsymbol{x}} = \left\{ \frac{\partial [\boldsymbol{A}(\boldsymbol{x}, \boldsymbol{\upmu}(t))]}{\partial \boldsymbol{x}} + \frac{\partial [\boldsymbol{B}(\boldsymbol{x}, \boldsymbol{\upmu}(t))]}{\partial \boldsymbol{x}} \right\}^{T} \boldsymbol{\upmu}(t+\varepsilon) - \\ &- \frac{1}{2} \frac{\partial [\boldsymbol{x}^{T}(t) \boldsymbol{Q}(\boldsymbol{x}, \boldsymbol{\upmu}(t))]}{\partial \boldsymbol{x}(t)} = \left\{ (\boldsymbol{x}^{T} \odot \boldsymbol{E}\_{\boldsymbol{n}})\_{\boldsymbol{n} \times \boldsymbol{n}^{2}} \left[ \frac{\partial \boldsymbol{A}(\boldsymbol{x}, \boldsymbol{\upmu})}{\partial \boldsymbol{x}} \right]\_{\boldsymbol{n}^{2} \times \boldsymbol{n}} + \boldsymbol{A}(\boldsymbol{x}, \boldsymbol{\upmu})\_{\boldsymbol{n} \times \boldsymbol{n}} + (\boldsymbol{u}^{T} \odot \boldsymbol{E}\_{\boldsymbol{n}})\_{\boldsymbol{n} \times \boldsymbol{n} \pi} \left[ \frac{\partial \boldsymbol{B}(\boldsymbol{x}, \boldsymbol{\upmu})}{\partial \boldsymbol{x}} \right]\_{\boldsymbol{n} \boldsymbol{n} \times \boldsymbol{n}} \right\}^{T} \boldsymbol{\upmu}(t+\varepsilon) - \\ &- \frac{1}{2} \boldsymbol{Q}(\boldsymbol{x}, \boldsymbol{\upmu})\_{\boldsymbol{n} \times \boldsymbol{n}} \boldsymbol{x} - \frac{1}{2} \left[ \frac{\partial \boldsymbol{Q}(\boldsymbol{x}, \boldsymbol{\upmu})}{\partial \boldsymbol{x}} \right]\_{\boldsymbol{n} \times \boldsymbol{n}^{2}}^{T} (\boldsymbol{x} \odot \boldsymbol{E}\_{\boldsymbol{n}})\_{\boldsymbol{n}$$

where ⊗ stands for Kronecker matrix product and *En* is a *n* × *n* unity matrix. From (6) we obtain the next expression for the control

$$
\mu(t) = R^{-1} B^T(\mathbf{x}, \boldsymbol{\varepsilon}) \psi(t + \boldsymbol{\varepsilon}).\tag{7}
$$

Using the representation *ψ*(*t*) = −*Pε*(*x*, *t*)*x*(*t*), we obtain the following expression instead of (7):

$$\mathbf{u}(\mathbf{x},t,\varepsilon) = -\left\{\mathbf{R} + \mathbf{B}^T(\mathbf{x},\varepsilon)\mathbf{P}(\mathbf{x},t+\varepsilon)\mathbf{B}(\mathbf{x},\varepsilon)\right\}^{-1}\mathbf{B}^T(\mathbf{x},\varepsilon)\mathbf{P}(\mathbf{x},t+\varepsilon)\mathbf{A}(\mathbf{x},\varepsilon)\mathbf{x}(t) = \mathbf{K}(\mathbf{x},\varepsilon,t+\varepsilon)\mathbf{x}(t),\tag{8}$$

where *P*(*x*, *t*,*ε*) must satisfy the discrete matrix Riccati-type equation with the zero-control matrix and as a result the missing quadratic nonlinearity part

$$\begin{split} & -P(\mathbf{x},t,\varepsilon) + \left\{ \left( \mathbf{x}^{T} \odot E\_{\mathrm{n}} \right)\_{\mathrm{n}\times\mathsf{R}^{2}} \left[ \frac{\partial A(\mathbf{x},t)}{\partial \mathbf{z}} \right]\_{\mathrm{n}^{2}\times\mathsf{R}} + A(\mathbf{x},\varepsilon)\_{\mathrm{n}\times\mathsf{R}} + \left( \mathbf{x}^{T}K(\mathbf{x},\varepsilon)^{T} \odot E\_{\mathrm{n}} \right)\_{\mathrm{n}\times\mathsf{R}\mathrm{T}} \left[ \frac{\partial B(\mathbf{x},t)}{\partial \mathbf{z}} \right]\_{\mathrm{n}\times\mathsf{R}\mathrm{T}} \right\}^{T} \times \\ & \times P(\mathbf{x},t+\varepsilon,\varepsilon) \times \\ & \times \left\{ E - B(\mathbf{x},\varepsilon) \left[ R + B^{T}(\mathbf{x},\varepsilon) P(\mathbf{x},t+\varepsilon,\varepsilon) B(\mathbf{x},\varepsilon) \right]^{-1} B^{T}(\mathbf{x},\varepsilon) P(\mathbf{x},t+\varepsilon,\varepsilon) \right\} A(\mathbf{x},\varepsilon) + \\ & + Q(\mathbf{x},\varepsilon)\_{\mathrm{n}\times\mathsf{R}} + \frac{1}{2} \left[ \frac{\partial Q(\mathbf{x},t)}{\partial \mathbf{z}} \right]\_{\mathrm{n}\times\mathsf{R}^{2}}^{T} \left( \mathbf{x} \odot E\_{\mathrm{n}} \right)\_{\mathrm{n}^{2}\times\mathsf{R}} = \Phi(P,t,\varepsilon) = 0, \end{split}$$

with the initial condition at the end of the interval *P*(*x*, *T*,*ε*) = *F*. The latter can be represented as

$$\begin{array}{l} \Phi(P,t,\varepsilon) = -P(\mathbf{x},t,\varepsilon) + A^T(\mathbf{x},\varepsilon)\_{n\times n}P(\mathbf{x},t+\varepsilon,\varepsilon)A(\mathbf{x},\varepsilon) - \\ -A^T(\mathbf{x},\varepsilon)\_{n\times n}P(\mathbf{x},t+\varepsilon,\varepsilon)B(\mathbf{x},\varepsilon)\left[R + B^T(\mathbf{x},\varepsilon)P(\mathbf{x},t+\varepsilon,\varepsilon)B(\mathbf{x},\varepsilon)\right]^{-1} \times \\ \times B^T(\mathbf{x},\varepsilon)P(\mathbf{x},t+\varepsilon,\varepsilon)A(\mathbf{x},\varepsilon) + Q(\mathbf{x},\varepsilon)\_{n\times n} + \varepsilon\Omega(P(\mathbf{x},t+\varepsilon),\mathbf{x},t+\varepsilon) = \\ = \Psi(P,t,\varepsilon) + \varepsilon\Omega(P(\mathbf{x},t+\varepsilon),\mathbf{x},t+\varepsilon) = 0, \ P(\mathbf{x},T,\varepsilon) = F,\end{array} \tag{9}$$

where Ψ(*P*, *t*,*ε*) is the usual discrete matrix Riccati operator and *<sup>K</sup>*(*x*,*ε*) = <sup>−</sup><sup>3</sup> *R* + *BT*(*x*,*ε*)*P*(*x*, *t* + *ε*,*ε*)*B*(*x*,*ε*) 4−<sup>1</sup> *BT*(*x*,*ε*)*P*(*x*, *t* + *ε*,*ε*)*A*(*x*,*ε*). The difference from the regular discrete difference Riccati equation on a finite interval

lies in an additional, but cumbersome, term on the right-hand side with the function

Ω(*P*(*x*, *t* + *ε*), *x*, *t* + *ε*) = = ⎧ ⎪⎪⎪⎪⎨ ⎪⎪⎪⎪⎩ (*x<sup>T</sup>* <sup>⊗</sup> *En*)*n*×*n*<sup>2</sup> \$ *<sup>∂</sup>A*1(*x*) *∂x* % *<sup>n</sup>*2×*<sup>n</sup>* − −(*x<sup>T</sup>* 3 *R* + *BT*(*x*,*ε*)*P*(*x*, *t* + *ε*,*ε*)*B*(*x*,*ε*) 4−<sup>1</sup> *BT*(*x*,*ε*)*P*(*x*, *t* + *ε*,*ε*)*A*(*x*,*ε*) !*T* ⊗ *En*) *<sup>n</sup>*×*nr*<sup>×</sup> × \$ *<sup>∂</sup>B*1(*x*) *∂x* % *nr*×*n* ⎫ ⎪⎪⎪⎪⎬ ⎪⎪⎪⎪⎭ *T* × ×*P*(*x*, *t* + *ε*,*ε*) 0 *<sup>E</sup>* <sup>−</sup> *<sup>B</sup>*(*x*,*ε*)[*<sup>R</sup>* <sup>+</sup> *<sup>B</sup>T*(*x*,*ε*)*P*(*x*, *<sup>t</sup>* <sup>+</sup> *<sup>ε</sup>*,*ε*)*B*(*x*,*ε*)]−<sup>1</sup> *BT*(*x*,*ε*)*P*(*x*, *t* + *ε*,*ε*) 1 *A*(*x*,*ε*)+ +<sup>1</sup> 2 \$ *<sup>∂</sup>*(*Q*1(*x*)+*εQ*2(*x*)) *∂x* %*T <sup>n</sup>*×*n*<sup>2</sup> (*<sup>x</sup>* <sup>⊗</sup> *En*)*n*2×*n*.

Let's denote the matrix of the closed-loop (*cl*) system for (1) as *Acl*(*x*,*ε*)=(*A*(*x*,*ε*) − *<sup>B</sup>*(*x*,*ε*)*R*−1*BT*(*x*,*ε*)*K*(*x*,*ε*))*x*.

Taking into account the dependency of matrices *Q*(*x*,*ε*), *A*(*x*,*ε*), *B*(*x*,*ε*) on the small parameter a uniform asymptotic approximation of the second order for the solution of the singularly perturbed initial problem (9) is constructed using the boundary layer functions method (BLFM) [14,15]. The solution is found in the reverse time in the following form:

$$P(\mathbf{x}, t, \varepsilon) = \overline{P}\_2(\mathbf{x}, t, \varepsilon) + \Pi\_2 P(\mathbf{x}, \tau, \varepsilon), \ \tau = (t - T)/\varepsilon = -1, -2, \dots, -N, \ N = T/\varepsilon\_r$$

where *P*2(*x*, *t*,*ε*) = *P*<sup>0</sup> + *εP*1(*x*, *t*) + *ε*2*P*2(*x*, *t*),

Π2*P*(*x*, *τ*,*ε*) = Π0*P*(*x*, *τ*) + *ε*Π1*P*(*x*, *τ*)) + *ε*2Π2*P*(*x*, *τ*)), are the second-order partial sums of the regular and boundary power series by parameter *ε*, respectively.

Now the operator Φ is transformed in the following way:

Φ(*P*, *t*,*ε*) = Φ *P*(*x*, *t*,*ε*), *t*,*ε* + Φ *P*(*x*, *τε* + *T*,*ε*) + Π*P*(*x*, *τ*,*ε*), *τ*,*ε* <sup>−</sup> <sup>Φ</sup> *P*(*x*, *τε* + *T*,*ε*), *τε* + *T*,*ε* = = Φ + ΠΦ,

> where Φ(*x*, *t*,*ε*) = Φ(*P*(*x*, *t*,*ε*), *t*,*ε*) = Φ0(*P*0, *t*) + *ε*Φ1(*P*0, *t*) + ... , ΠΦ(*x*, *τε* + *T*,*ε*) = Φ *P*(*x*, *τε* + *T*,*ε*) + Π*P*(*x*, *τ*,*ε*), *τε* + *T*,*ε* <sup>−</sup> <sup>Φ</sup> *P*(*x*, *τε* + *T*,*ε*), *τε* + *T*,*ε* . Substituting the representations for *P*2(*x*, *t*,*ε*), Π2*P*(*x*, *τ*,*ε*) into the equation and the initial condition in (9) and equating the terms with the same powers of the small parameter, we obtain a discrete algebraic Riccati equation (DARE) for the zero term of the regular series *P*<sup>0</sup>

$$A\_0 \, ^T \overline{P}\_0 A\_0 - \overline{P}\_0 - A\_0 \, ^T \overline{P}\_0 B\_0 (R\_0 + B\_0 \, ^T \overline{P}\_0 B\_0) \, ^{-1} B\_0 \, ^T \overline{P}\_0 A\_0 + Q\_0 = 0,\tag{10}$$

and for the first term of the regular series *P*1(*x*), we obtain the matrix discrete algebraic Lyapunov equation

$$A\_{cl}^{0}{}^{T}\overline{P}\_{1}(\mathbf{x})A\_{cl}^{0} - \overline{P}\_{1}(\mathbf{x}) = -G\_{1}(\mathbf{x}),\tag{11}$$

where *A*<sup>0</sup> *cl* = *A*<sup>0</sup> + *B*0*K*<sup>0</sup> = *A*<sup>0</sup> − *B*0[*R*<sup>0</sup> + *B*<sup>0</sup> *TP*0*B*0] −1 *B*0 *TP*0*A*0− is the matrix of the linear closed-loop system (*cl*) and *G*1(*x*) = *D*0(*x*) + *D*1(*x*) − *D*2(*x*) − *D*3(*x*) + *D*4(*x*),

*D*0(*x*) = \$ *<sup>∂</sup>A*1(*x*) *∂x* %*T* (*x* ⊗ *En*) − \$ *∂B*1(*x*) *∂x* %*T* ( '*R* + *B*<sup>0</sup> *TP*0*B*0)−1*B*<sup>0</sup> *TP*0*A*<sup>0</sup> ( *x* ⊗ *En* 2 *P*0[*A*<sup>0</sup> − *B*0[ *R* + *B*<sup>0</sup> *TP*0*B*0] <sup>−</sup>1*B*<sup>0</sup> *TP*0*A*0]+ +*Q*1(*x*) + <sup>1</sup> 2 \$ *∂Q*1(*x*) *∂x* %*T n*×*n*<sup>2</sup> (*<sup>x</sup>* <sup>⊗</sup> *En*)*n*<sup>2</sup>×*n*, *D*1(*x*) = *A*<sup>0</sup> *TP*0*A*<sup>1</sup> + *A*<sup>1</sup> *<sup>T</sup>*(*x*)*P*0*A*0, *D*2(*x*) = *A*<sup>0</sup> *TP*0*B*<sup>1</sup> ' *R* + *B*<sup>0</sup> *TP*0*B*0] <sup>−</sup>1*B*<sup>0</sup> *TP*0*A*<sup>0</sup> + *A*<sup>0</sup> *TP*0*B*0[ *R* + *B*<sup>0</sup> *TP*0*B*0] <sup>−</sup>1*B*<sup>1</sup> *TP*0*A*0, *D*3(*x*) = *A*<sup>0</sup> *TP*0*B*<sup>0</sup> ' *R* + *B*<sup>0</sup> *TP*0*B*0] <sup>−</sup>1*B*<sup>0</sup> *TP*0*A*<sup>1</sup> + *A*<sup>1</sup> *<sup>T</sup>*(*x*)*P*0*B*0[ *R* + *B*<sup>0</sup> *TP*0*B*0] <sup>−</sup>1*B*<sup>0</sup> *TP*0*A*0, *D*4(*x*) = *A*<sup>0</sup> *TP*0*B*<sup>0</sup> ' *R* + *B*<sup>0</sup> *TP*0*B*0] <sup>−</sup>1*B*<sup>0</sup> *TP*0*B*1[ *R* + *B*<sup>0</sup> *TP*0*B*0] <sup>−</sup>1*B*<sup>0</sup> *TP*0*A*0+ +*A*<sup>0</sup> *TP*0*B*<sup>0</sup> ' *R* + *B*<sup>0</sup> *TP*0*B*0] <sup>−</sup>1*B*<sup>1</sup> *TP*0*B*0[ *R* + *B*<sup>0</sup> *TP*0*B*0] <sup>−</sup>1*B*<sup>0</sup> *TP*0*A*0.

For the second term of the regular series *P*2(*x*), the matrix discrete algebraic Lyapunov equation has the same form as Equation (11) for *P*1(*x*),

$$A\_{cl}^{0}{}^{T}\vec{P}\_{2}(\mathfrak{x})A\_{cl}^{0} - \vec{P}\_{2}(\mathfrak{x}) = -G\_{2}(\mathfrak{x})\_{\prime\prime}$$

where *G*2(*x*) has a similar structure as *G*1(*x*) but it is a more complex function of the found expansion terms and matrices of the system, and we omit its representation here.

For the zero term of the boundary series Π0*P*(*τ*), we obtain the difference initial problem

$$\begin{cases} \Pi\_{0}P(\tau) = -\overline{P}\_{0} + A\_{0}{}^{T}(\overline{P}\_{0} + \Pi\_{0}P(\tau+1))A\_{0} - A\_{0}{}^{T}(\overline{P}\_{0} + \Pi\_{0}P(\tau+1))B\_{0} \times \\ \times \left[\overline{R} + B\_{0}{}^{T}(\overline{P}\_{0} + \Pi\_{0}P(\tau+1))B\_{0}\right]^{-1}B\_{0}{}^{T}(\overline{P}\_{0} + \Pi\_{0}P(\tau+1))A\_{0} + Q\_{0\prime} \\ \Pi\_{0}P(0) = F - \overline{P}\_{0\prime} \end{cases} \tag{12}$$

and for the first term of the boundary series Π1*P*(*x*, *τ*), we have the following discrete problem:

$$\begin{array}{l} \Pi\_{1}P(\mathbf{x},\tau) = \widetilde{\zeta}^{1}(\tau)\Pi\_{1}P(\tau+1)\widetilde{\zeta}(\tau) + \Pi\_{1}G(\mathbf{x},\tau),\\ \Pi\_{1}P(0) = -\widetilde{P}\_{1}(\mathbf{x}),\ \tau = -1, -2, \dots, \end{array} \tag{13}$$

where *ζ*(*τ*) = *A*<sup>0</sup> − *B*0[*R* + *B*<sup>0</sup> *<sup>T</sup>*(*P*<sup>0</sup> + Π0*P*(*τ* + 1))*B*0] −1 *B*0 *<sup>T</sup>*(*P*<sup>0</sup> + Π0*P*(*τ* + 1))*A*0, and Π1*G*(*x*, *τ*) has a complex structure

Π1*G*(*x*, *τ*) = *A*<sup>0</sup> *<sup>T</sup>*(*P*<sup>0</sup> + Π0*P*(*τ* + 1))*B*0(*R*<sup>0</sup> + *B*<sup>0</sup> *<sup>T</sup>*(*P*<sup>0</sup> + Π0*P*(*τ* + 1))*B*0) −1 *B*0 *TP*1(*x*)*B*0× × *R*<sup>0</sup> + *B*<sup>0</sup> *T P*<sup>0</sup> + Π0*P*(*τ* + 1) *B*0 −<sup>1</sup> *B*0 *T P*<sup>0</sup> + Π0*P*(*τ* + 1) *A*0− −*A*<sup>0</sup> *<sup>T</sup>*(*P*<sup>0</sup> + Π0*P*(*τ* + 1))*B*0(*R*<sup>0</sup> + *B*<sup>0</sup> *<sup>T</sup>*(*P*<sup>0</sup> + Π0*P*(*τ* + 1))*B*0) −1 *B*0 *TP*1(*x*)*A*<sup>0</sup> +*A*<sup>0</sup> *TP*1(*x*)*A*<sup>0</sup> − *<sup>B</sup>*0(*R*<sup>0</sup> + *<sup>B</sup>*<sup>0</sup> *<sup>T</sup>*(*P*<sup>0</sup> + Π0*P*(*τ* + 1))*B*0) −1 *B*0 *<sup>T</sup>*(*P*<sup>0</sup> + Π0*P*(*τ* + 1))*A*0)+ −*P*1(*x*) + *A*<sup>0</sup> *T P*<sup>0</sup> + Π0*P*(*τ* + 1) *A*1(*x*)+ +*A*<sup>0</sup> *T P*<sup>0</sup> + Π0*P*(*τ* + 1) (−*B*<sup>0</sup> *R*<sup>0</sup> + *B*<sup>0</sup> *T P*<sup>0</sup> + Π0*P*(*τ* + 1) *B*0 −<sup>1</sup> × × *B*1 *<sup>T</sup>*(*x*) *P*<sup>0</sup> + Π0*P*(*τ* + 1) *A*<sup>0</sup> + *B*<sup>0</sup> *T P*<sup>0</sup> + Π0*P*(*τ* + 1) *A*1(*x*) + +*B*<sup>0</sup> *R*<sup>0</sup> + *B*<sup>0</sup> *T P*<sup>0</sup> + Π0*P*(*τ* + 1) *B*0 −<sup>1</sup> × × *B*1 *<sup>T</sup>*(*x*) *P*<sup>0</sup> + Π0*P*(*τ* + 1) *B*<sup>0</sup> + *B*<sup>0</sup> *T P*<sup>0</sup> + Π0*P*(*τ* + 1) *B*1(*x*) × × *R*<sup>0</sup> + *B*<sup>0</sup> *T P*<sup>0</sup> + Π0*P*(*τ* + 1) *B*0 −<sup>1</sup> *B*0 *T P*<sup>0</sup> + Π0*P*(*τ* + 1) *A*0− −*B*1(*x*) *R*<sup>0</sup> + *B*<sup>0</sup> *T P*<sup>0</sup> + Π0*P*(*τ* + 1) *B*0 −<sup>1</sup> *B*0 *T P*<sup>0</sup> + Π0*P*(*τ* + 1) *A*0)+ +(\$ *<sup>∂</sup>A*1(*x*) *∂x* %*T* (*x* ⊗ *En*) + *A*<sup>1</sup> *<sup>T</sup>*(*x*)− − \$ *<sup>∂</sup>B*1(*x*) *∂x* %*<sup>T</sup> R*<sup>0</sup> + *B*<sup>0</sup> *T P*<sup>0</sup> + Π0*P*(*τ* + 1) *B*0 −<sup>1</sup> *B*0 *T P*<sup>0</sup> + Π0*P*(*τ* + 1) *A*<sup>0</sup> ! *x* ⊗ *En* ! )× × *P*<sup>0</sup> + Π0*P*(*τ* + 1) *A*<sup>0</sup> − *B*<sup>0</sup> *R*<sup>0</sup> + *B*<sup>0</sup> *T P*<sup>0</sup> + Π0*P*(*τ* + 1) *B*0 −<sup>1</sup> *B*0 *T P*<sup>0</sup> + Π0*P*(*τ* + 1) *A*<sup>0</sup> ! .

For Π2*P*(*x*, *τ*), we obtain a discrete initial problem similar to (12), but Π2*G*(*x*, *τ*) has an even more cumbersome structure than Π1*G*(*x*, *τ*) and we omit this representation here. Let us introduce the conditions:

*I. Coefficients of matrices A*1(*x*), *B*1(*x*), *Q*1(*x*), *Q*2(*x*) *are continuous and bounded functions on* Х*, g*(0,*ε*) ≡ 0 *(g*(*x*,*ε*) = *Acl*(*x*(*t*),*ε*)*x*(*t*)*) and the parameterεbelongs to a bounded interval* (0,*ε*0]*, the solution to problems (1) and (2) exists and is bounded for all admissible x*0,*ε*, *t*; 1

*II. The triple of matrices* (*A*0, *B*0, *Q*<sup>0</sup> <sup>2</sup> )*is controllable and observable.*

By definition, a system is said to be controllable, if it is possible to transfer the system state from any initial state to any desired state within a finite interval of time. A system is said to be observable if every state can be completely identified by measurements of the outputs at a finite time interval.

The following is true.

**Theorem 1.** *Under conditions I–II, there is a sufficiently small value of ε*<sup>0</sup> > 0*, such that for all x* ∈ *X*, *t* ∈ [0, *T*], *ε* ∈ (0,*ε*0], *τ* = −1, −2, . . . , −*N*,... *the following statements hold:*


$$\left\|\left\|P(\mathbf{x},t,\varepsilon)-(\overline{P}\_{2}(\mathbf{x},t,\varepsilon)+\Pi\_{2}P(\mathbf{x},\tau,\varepsilon))\right\|\right\|=O(\varepsilon^{3})\\\begin{cases} \varepsilon=O(\varepsilon^{3})\\ \mathbf{x}\in X, \ t\in[0,T], \ \varepsilon\in(0,\varepsilon\_{0}], \ \tau=-1,-2,\ldots,-N,\ldots \end{cases} \tag{14}$$

**Proof of Theorem 1.** The statements of the theorem for each *x* ∈ *X* follow from the corresponding proof schemes of the statements in paragraph 2 of paper [18], where the similar linear-quadratic problem is considered, but with time-dependent matrix coefficients. The statements here generalize the similar results presented in [18] since in comparison with [18] the initial problem (9) is additionally regularly perturbed with continuously differentiable perturbation *ε*Ω on the right-hand side of the matrix discrete Riccati equation. Therefore, here we present only the new components of the proof.

As the associated system for the singularly perturbed problem (9) for *ε* = 0, Ω ≡ 0 coincides with the analogous one in [18], where it is shown that the positive definite root *P*<sup>0</sup> of the limiting algebraic discrete matrix Riccati equation is an asymptotically stable equilibrium point of system (9) for *τ* → −∞ and matrix *F* belongs to the domain of influence of this root. Further, in [18], the form of the main functional matrix Γ of the associated system to problem (9) is established, which is calculated by setting *ε* = 0 in the right-hand side and using the Kronecker product of matrices can be represented in the form Γ = *A*<sup>0</sup> *cl <sup>T</sup>* ⊗ *<sup>A</sup>*<sup>0</sup> *cl*. From the properties of the spectrum of the constant matrix *<sup>A</sup>*<sup>0</sup> *cl* (the matrix of the corresponding closed-loop system) under condition *II*, it follows that its eigenvalues *λi*, *i* = 1, 2, ... , *n* are inside the unit circle. Taking into account the properties of the spectrum of the Kronecker product of matrices, we find that the eigenvalues of matrix Γ are representable in the form *λiλj*, *i*, *j* = 1, 2, ... , *n*. Thus, it is established that the spectrum of the main functional matrix of the associated system has eigenvalues inside the unit circle.

Because the presence of a disturbance *ε*Ω will not fundamentally change the form of inhomogeneities in the corresponding Riccati and Lyapunov equations, then, as in Theorem 2.1 from [18], it can be proved that there exist such *α* > 0, *β* > 0, that the corresponding estimates for the boundary functions hold

$$\|\Pi\_i P(\tau)\| \le \kappa \exp(\beta \tau), \ \tau \le 0, \ i = 0, 1, 2. \tag{15}$$

Despite the presence of regular perturbations *ε*Ω in Equation (9) the next estimate for the residual term of the asymptotics *<sup>η</sup>*(*x*, *<sup>t</sup>*, *<sup>τ</sup>*,*ε*) = *<sup>P</sup>*(*x*, *<sup>t</sup>*,*ε*) <sup>−</sup> <sup>2</sup> ∑ *i*=0 *εi* (*Pi*(*x*, *t*) + Π*iP*(*τ*)), *τ* = *<sup>t</sup>*−*<sup>T</sup> <sup>ε</sup>* can be obtained using the scheme from [18],

$$\|\|\eta(\mathbf{x}, \mathbf{t}, \tau, \varepsilon)\|\| = O(\varepsilon^3), \tag{16}$$

for all *x* ∈ *X*, *ε* ∈ (0,*ε*0], *t* ∈ [0, *T*], *τ* = −1, −2, . . . , −*N*.

Moreover, using the smoothness of all functions included in Φ, Ω and assuming that the norm *P*(*x*, *t*,*ε*) is uniformly bounded for all *x* ∈ *X*, *t* ∈ [0, *T*], *ε* ∈ (0,*ε*0], *τ* = −1, −2, ... , −*N*, ..., here we can follow the proof as in Lemma 6.1 and Theorem 6.1 given in [14].

This, in turn, will lead to the fact that when choosing a sufficiently small *ε*<sup>0</sup> > 0, one can obtain the existence, as well as the uniqueness of the solution in the problem (9). The last statements follow from the application of the principle of contractive mappings to the equation for the residual term of asymptotics *η*(*x*, *t*, *τ*,*ε*).

#### *2.2. Symmetrization*

Since the matrix *G*1(*x*) may not be symmetric, we introduce the transposed equation, and by adding the equation *A*<sup>0</sup> *cl TP*1(*x*)*A*<sup>0</sup> *cl* − *<sup>P</sup>*1(*x*) = *<sup>G</sup>*1(*x*)*A*<sup>0</sup> *cl TP*1(*x*)*A*<sup>0</sup> *cl* − *P*1(*x*) = *G*1(*x*) to *A*<sup>0</sup> *cl TP*<sup>1</sup> *<sup>T</sup>*(*x*)*A*<sup>0</sup> *cl* − *P*<sup>1</sup> *<sup>T</sup>*(*x*) = *G*<sup>1</sup> *<sup>T</sup>*(*x*), we get

$$A\_{cl}^{0,T}(\overline{P}\_1(\mathbf{x}) + \overline{P}\_1^{,T}(\mathbf{x}))A\_{cl}^0 - (\overline{P}\_1(\mathbf{x}) + \overline{P}\_1^{,T}(\mathbf{x})) = G\_1(\mathbf{x}) + G\_1^{,T}(\mathbf{x})\_{\perp}$$

Further, in this paper the next equations with "averaged" right-hand sides will be used for the regular and the boundary layer terms of the first and second order (the terms of the power series *ε* and *ε*2). For example, instead of Equation (11) the following equation is introduced for a symmetric matrix ) *P*1(*x*)

$$A\_{cl}^{0}{}^{T}\dot{\overline{P}}\_{1}(\mathbf{x})A\_{cl}^{0} - \dot{\overline{P}}\_{1}(\mathbf{x}) = \widetilde{G}\_{1}(\mathbf{x}),\tag{17}$$

where *<sup>G</sup>*)1(*x*) = <sup>1</sup> 2 *G*1(*x*) + *G*<sup>1</sup> *<sup>T</sup>*(*x*) , *<sup>P</sup>*)1(*x*, *<sup>t</sup>*) = <sup>1</sup> <sup>2</sup> () *P*1(*x*) + ) *P*1 *<sup>T</sup>*(*x*)). A similar operation is performed for the second approximation term.

In addition, Π) <sup>1</sup>*P*(*x*, *τ*) is found in the form

$$\begin{array}{l} \stackrel{\cdot}{\Pi}\_{1}P(\mathbf{x},\mathsf{\tau}) = \zeta^{T}(\mathsf{x},\mathsf{\tau})\dot{\Pi}\_{1}P(\mathsf{\tau}+1)\zeta(\mathsf{x},\mathsf{\tau}) + \dot{\Pi}\_{1}G(\mathsf{x},\mathsf{\tau}),\\ \stackrel{\cdot}{\Pi}\_{1}P(0) = -\overline{P}\_{1}(\mathsf{x}),\ \mathsf{\tau} = -1, -2, \dots, \end{array}$$

where <sup>Π</sup>) <sup>1</sup>*G*(*x*) = <sup>1</sup> 2 Π1*G*(*x*) + Π1*GT*(*x*) , <sup>Π</sup>) <sup>1</sup>*P*(*x*, *<sup>τ</sup>*) = <sup>1</sup> <sup>2</sup> (Π1*P*(*x*, *<sup>τ</sup>*) + <sup>Π</sup>1*PT*(*x*, *<sup>τ</sup>*)).

A similar equation is obtained for Π) <sup>2</sup>*P*(*x*, *τ*), where the inhomogeneity in the righthand side equals to <sup>Π</sup>) <sup>2</sup>*G*(*x*) = <sup>1</sup> 2 Π2*G*(*x*) + Π2*GT*(*x*) .

#### **3. Discrete One-Point Padé Regulator**

The asymptotic analysis for small parameter values can lead to an acceptable quality approximation of the exact solution, but with an increase of the small parameter value the asymptotic representations can strongly deviate from the exact solutions and their use in numerical analysis for larger values of the parameter is limited and at best they can serve only to restore the qualitative nature of the solution behavior. *PA* often increases the interval of parameter variation for which it can provide the approximation of the exact solution and restore its qualitative picture in comparison with the asymptotics. Thus, *PA* demonstrates better extrapolation properties [23].

In general, a particular system of algebraic equations, which, generally speaking, is in a certain way selected from some redefined system, is solved to find the *PA*.

Here a one-point Padé regulator of an order [1/2] is constructed, which contains two asymptotic approximations: the first order uniform asymptotic approximation in the "numerator" of the *PA*—for reproducing the boundary layer in the general construction, and the second-order approximation of some regular series in the "denominator", i.e., the proposed construction has the following form:

$$\begin{array}{l} \text{PA}\_{[1/2]}(\mathbf{x}, t, \tau, \varepsilon) = \left( \mathbf{M}\_{0}(\mathbf{x}) + \varepsilon \mathbf{M}\_{1}(\mathbf{x}) + \Gamma \mathbf{M}\_{0}(\mathbf{x}, \tau) + \varepsilon \Gamma \mathbf{M}\_{1}(\mathbf{x}, \tau) \right) \times \\ \times \left( E + \varepsilon \mathbf{N}\_{1}(\mathbf{x}) + \varepsilon^{2} \mathbf{N}\_{2}(\mathbf{x}) \right)^{-1} . \end{array} \tag{18}$$

Note that the form of the "denominator" in (18) is less complex than the "numerator" which makes it easier to overcome the "denominator" zeros problem, which is the main reason for the quality decline of the approximations of the exact solution using *PA*.

So, we do not introduce the boundary functions in the "denominator" of *PA*. Decomposing the matrix in (18) in a series of integer powers of parameter *ε* and equating the terms with the same powers of the parameter *ε* of the resulting decomposition and the corresponding terms of the expansions *φ* and ΠΦ, separately for the terms dependent on *t* and *τ*, we get an inhomogeneous linear system of six equations with matrix coefficients

depending on *x* ∈ *X*, *t* ∈ [0, *T*], *ε* ∈ (0,*ε*0], *τ* = −1, −2, ... , −*N* for the six unknown matrices, the coefficients of *PA*.

$$\begin{array}{l} \mathsf{M}\_{0}(\mathbf{x}) = \mathsf{P}\_{0} \\ \mathsf{I}\,\mathrm{IM}\_{0}(\mathbf{x},\tau) = \Pi\_{0}\mathrm{P}(\tau) \\ \mathsf{M}\_{1}(\mathbf{x}) - \mathsf{P}\_{0}\mathrm{N}\_{1}(\mathbf{x}) = \mathsf{P}\_{1}(\mathbf{x}) \\ \mathsf{I}\,\mathrm{IM}\_{1}(\mathbf{x},\tau) - \Pi\_{0}\mathrm{P}(\tau)\mathrm{N}\_{1}(\mathbf{x}) = 0 \\ \mathsf{P}\_{1}(\mathbf{x})\mathrm{N}\_{1}(\mathbf{x}) + \mathsf{P}\_{0}\mathrm{N}\_{2}(\mathbf{x}) = -\mathsf{P}\_{2}(\mathbf{x}) \\ \Pi\_{1}\mathrm{P}(\mathbf{x},\tau)\mathrm{N}\_{1}(\mathbf{x}) + \Pi\_{0}\mathrm{P}(\tau)\mathrm{N}\_{2}(\mathbf{x}) = -\Pi\_{1}\mathrm{P}(\mathbf{x},\tau). \end{array} \tag{19}$$

The first matrices M0(*x*), ΠM0(*x*, *τ*) are immediately determined from the first two equations, and for the remaining four matrices the following linear system is obtained:

$$
\begin{pmatrix} E\_{\pi} & 0 & -\overline{\mathsf{P}}\_{0} & 0 \\ 0 & E\_{\pi} & -\Pi\_{0}P(\mathsf{r}) & 0 \\ 0 & 0 & \overline{\mathsf{P}}\_{1}(\mathsf{x}) & \overline{\mathsf{P}}\_{0} \\ 0 & 0 & \Pi\_{1}P(\mathsf{x},\mathsf{r}) & \Pi\_{0}P(\mathsf{r}) \end{pmatrix} \begin{pmatrix} \mathcal{M}\_{1}(\mathsf{x}) \\ \Pi\mathcal{M}\_{1}(\mathsf{x},\mathsf{r}) \\ \mathcal{N}\_{1}(\mathsf{x}) \\ \mathcal{N}\_{2}(\mathsf{x}) \end{pmatrix} = \begin{pmatrix} \overline{\mathsf{P}}\_{1}(\mathsf{x}) \\ 0 \\ -\overline{\mathsf{P}}\_{2}(\mathsf{x}) \\ -\Pi\mathcal{P}(\mathsf{x},\mathsf{r}) \end{pmatrix},\tag{20}
$$

where *En*—is an identity matrix of the size *n* × *n*. Next by denoting the matrix of a linear system (20) by *Y* = *Y*<sup>11</sup> *Y*<sup>12</sup> *<sup>Y</sup>*<sup>21</sup> *<sup>Y</sup>*<sup>22</sup> , where the corresponding blocks are *Y*<sup>11</sup> = *E*2*n*, *Y*<sup>21</sup> = 0, *Y*<sup>12</sup> = <sup>−</sup>*P*<sup>0</sup> <sup>0</sup> −Π0*P*(*τ*) 0 , *Y*<sup>22</sup> = *P*1(*x*) *P*<sup>0</sup> Π1*P*(*x*, *τ*) Π0*P*(*τ*) and by taking into account that block *Y*<sup>11</sup> = *E*2*<sup>n</sup>* is a nondegenerate (2*n* × 2*n*) identity matrix, we find that for the existence of the matrix inverse *Y*−<sup>1</sup> = *<sup>E</sup>* <sup>−</sup>*Y*12*Y*22−<sup>1</sup> 0 Δ−<sup>1</sup> it is necessary and sufficient [24] that matrix <sup>Δ</sup> <sup>=</sup> *<sup>Y</sup>*<sup>22</sup> <sup>−</sup> *<sup>Y</sup>*21*Y*−<sup>1</sup> <sup>11</sup> *Y*<sup>12</sup> is nondegenerate and Δ = *Y*<sup>22</sup> since *Y*<sup>21</sup> ≡ 0. Now the solution of the last two equations in (20) is explicitly defined by

$$
\Gamma \begin{pmatrix} N\_1(\mathbf{x}) \\ N\_2(\mathbf{x}) \end{pmatrix} = \boldsymbol{\Delta}^{-1} \begin{pmatrix} -\mathsf{P}\_2(\mathbf{x}) \\ -\Pi\_2 P(\mathbf{x}, \mathbf{r}) \end{pmatrix}. \tag{21}
$$

As M1(*x*) and ΠM1(*x*, *τ*) are found from the first two equations in (20) and are the functions of *N*1(*x*), where *N*1(*x*), *N*2(*x*) are determined from the last two equations in (20) and take the form (21).

As some of the matrices in (18)–(20) that are the regular series terms in the asymptotic representation of *P*(*x*, *t*,*ε*) are positive definite, it is possible to make the other matrices sign-definite to guarantee the solvability of (20) and positive definiteness of *PA* on the entire interval, for example, by a special choice of criteria matrices *F*, *Q*<sup>0</sup> > 0, *Q*1(*x*) > 0, *Q*2(*x*) > 0, *x* ∈ *X*. In the latter case, we are dealing not with a problem of the criterion (2) minimization along the trajectories of (1), but with a synthesis construction problem using the SDRE algorithm.

Let us introduce a modified *PA*, defined as *PA*) [1/2](*x*, *t*, *τ*,*ε*), for which the system of equations is obtained from system (20) by replacing the elements of the first and secondorder terms of the approximation with their symmetric, "averaged" values.

$$
\begin{pmatrix} E\_{\pi} & 0 & -\overline{P}\_{0} & 0 \\ 0 & E\_{\pi} & -\Pi\_{0}P(\pi) & 0 \\ 0 & 0 & \widetilde{\mathcal{P}}\_{1}(\mathbf{x}) & \overline{P}\_{0} \\ 0 & 0 & \widetilde{\Pi}\_{1}P(\mathbf{x},\boldsymbol{\tau}) & \Pi\_{0}P(\boldsymbol{\tau}) \end{pmatrix} \begin{pmatrix} \mathcal{M}\_{1}(\mathbf{x}) \\ \Pi\mathcal{M}\_{1}(\mathbf{x},\boldsymbol{\tau}) \\ \mathcal{N}\_{1}(\mathbf{x}) \\ \mathcal{N}\_{2}(\mathbf{x}) \end{pmatrix} = \begin{pmatrix} \widetilde{\mathcal{P}}\_{1}(\mathbf{x}) \\ 0 \\ -\widetilde{\mathcal{P}}\_{2}(\mathbf{x}) \\ -\widetilde{\Pi}\_{2}P(\mathbf{x},\boldsymbol{\tau}) \end{pmatrix} \tag{22}
$$

By analogy with the study of system (20), the last two equations are firstly solved and the following matrices are formally introduced, Δ) = *Y*)<sup>22</sup> = <sup>6</sup> ) *P*1(*x*) *P*<sup>0</sup> Π) <sup>1</sup>*P*(*x*, *τ*) Π0*P*(*τ*) 7 and *Z* = Π0*P* − Π) <sup>1</sup>*P*(*x*, *τ*)) *P*1(*x*) −1 *P*0, *L* = ) *P*1(*x*) − *P*0Π0*P*(*τ*) −1 Π) <sup>1</sup>*P*(*x*, *τ*), which, in the case of their non-degeneracy allow us to present the solution for the "denominator" matrices in an explicit form *<sup>N</sup>*1(*x*) *N*2(*x*) = 6 −*L*−1) *P*2(*x*) + ) *P*1(*x*) −1 *<sup>P</sup>*0*Z*−1Π) <sup>2</sup>*P*(*x*, *<sup>τ</sup>*) *<sup>Z</sup>*−1[Π) <sup>1</sup>*P*(*x*, *<sup>τ</sup>*)) *P*1(*x*) <sup>−</sup>1) *<sup>P</sup>*2(*x*) <sup>−</sup> <sup>Π</sup>) <sup>2</sup>*P*(*x*, *<sup>τ</sup>*)] <sup>7</sup> , <sup>Δ</sup>)−<sup>1</sup> = 6 *<sup>K</sup>*−<sup>1</sup> −) *P*1(*x*) −1 *P*0*H*−<sup>1</sup> −*H*−1Π) <sup>1</sup>*P*(*x*, *<sup>τ</sup>*)) *P*1(*x*) <sup>−</sup><sup>1</sup> *H*−<sup>1</sup> 7 .

For consistency with the Kalman regulator (the linear-quadratic regulator named after R.E. Kalman who posed and solved the corresponding control problem for nonstationary linear systems in 1960), which leads to stabilization and also to the optimal trajectory according to the quadratic criterion in a closed-loop linear system on the semi-axis, the following symmetric gain matrix is introduced to ensure the symmetry of the gain matrix of the regulator based on the *PA* for *P*(*x*, *t*,*ε*):

$$K\_{\left[1/2\right]}(\mathbf{x},\mathbf{\tau},\varepsilon) = \frac{1}{2} (\check{P}A\_{\left[1/2\right]}(\mathbf{x},\mathbf{\tau},\varepsilon) + \check{P}A\_{\left[1/2\right]}{}^{T}(\mathbf{x},\mathbf{\tau},\varepsilon)),\tag{23}$$

which leads to a modified Padé regulator for the *PA* obtained from the system (22)

$$u(\mathbf{x},t,\varepsilon) = -\left\{\mathbf{R} + B^T(\mathbf{x},\varepsilon)\mathbf{K}\_{[1/2]}(\mathbf{x},\tau,\varepsilon)B(\mathbf{x},\varepsilon)\right\}^{-1}B^T(\mathbf{x},\varepsilon)\mathbf{K}\_{[1/2]}(\mathbf{x},\tau,\varepsilon)A(\mathbf{x},\varepsilon)\mathbf{x}(t). \tag{24}$$

**Remark 1.** *In numerical experiments, it becomes possible to modify the proposed PA structure (18) by the introduction of multipliers in front of the matrix coefficients in systems (19)–(20) and search for the values of these multipliers.*

The following conditions are additionally introduced

*III. Matrices F*, *Q*<sup>0</sup> > 0, *Q*1(*x*) > 0, *Q*2(*x*) > 0 *can be selected such that F* − *P*<sup>0</sup> > 0*, <sup>G</sup>*)1(*x*) > 0, *<sup>G</sup>*)2(*x*) > 0, *matrix* <sup>Δ</sup>)−<sup>1</sup> *exists and is uniformly bounded and* <sup>Π</sup>0*P*(*τ*) > <sup>0</sup> *for all x* ∈ *X*, *t* ∈ [0, *T*], *ε* ∈ (0,*ε*0], *τ* = −1, −2, . . . , 1 − *T*.

*IV. Matrix E* + *εN*1(*x*) + *ε*2*N*2(*x*) −<sup>1</sup> *exists, where*

$$\begin{cases} N\_1(\mathbf{x}) = -L^{-1}\overleftarrow{\mathcal{P}}\_2(\mathbf{x}) + \overleftarrow{\mathcal{P}}\_1(\mathbf{x})^{-1}\overleftarrow{\mathcal{P}}\_0\texttt{Z}^{-1}\overleftarrow{\mathcal{I}}\_2\mathrm{P}(\mathbf{x},\tau),\\ N\_2(\mathbf{x}) = Z^{-1}[\overleftarrow{\mathcal{I}}\_1\mathrm{P}(\mathbf{x},\tau)\overleftarrow{\mathcal{P}}\_1(\mathbf{x})^{-1}\overleftarrow{\mathcal{P}}\_2(\mathbf{x}) - \overleftarrow{\mathcal{I}}\_2\mathrm{P}(\mathbf{x},\tau)],\\ Z = \Pi\_0\mathrm{P}(\tau) - \overleftarrow{\mathcal{I}}\overleftarrow{\mathcal{I}}\_1\mathrm{P}(\mathbf{x},\tau)\overleftarrow{\mathcal{P}}\_1(\mathbf{x})^{-1}\overleftarrow{\mathcal{P}}\_0,\ L = \overleftarrow{\mathcal{P}}\_1(\mathbf{x}) - \overleftarrow{\mathcal{P}}\_0\Pi\_0\mathrm{P}(\tau)^{-1}\overleftarrow{\mathcal{I}}\_1\mathrm{P}(\mathbf{x},\tau) \end{cases}$$

for all *x* ∈ *X*, *t* ∈ [0, *T*], *ε* ∈ (0,*ε*0], *τ* = −1, −2, . . . , 1 − *T*.

For the construction of *PA* [1/2] and the solvability of the corresponding system for the coefficients of *PA* an asymptotic expansion of the second order is required, where some of the terms, in particular of the second order, can be found approximately.

The next statement takes place.

**Theorem 2.** *If conditions I–IV are satisfied, then there is a sufficiently small ε*<sup>0</sup> > 0*, such that for all x* ∈ *X*, *t* ∈ [0, *T*], *ε* ∈ (0,*ε*0], *τ* = −1, −2, ... , 1 − *T there is a unique solution of the matrix system of Equations (20) and the corresponding one-point matrix PA [1/2] of form (18) with a symmetric matrix K*[1/2](*x*, *τ*,*ε*) *(23) exists.*

**Proof of Theorem 2.** From condition II *P*<sup>0</sup> > 0 we get M0(*x*) > 0, then by condition III Π0*P*(*τ*) > 0 ∀*τ* and from here ΠM0(*x*, *τ*) > 0. Let us consider the discrete linear Lyapunov Equation (17) for ) *P*1(*x*). It is known [24] that if *G*)1(*x*) > 0 (condition III), the solution ) *P*1(*x*) of this equation is positive definite. From condition III, *G*)2(*x*) > 0 and there exist ) *P*2(*x*) > 0. It is easy to show that *ζ*(*x*, *τ*), Π) <sup>1</sup>*G*(*x*) are found from (13) and the corresponding matrices Π) <sup>1</sup>*P*(*x*, *τ*) are obtained as the solutions of difference Lyapunov equations with the initial condition Π) <sup>1</sup>*P*(0) = −*P*1(*x*), *τ* = −1, −2, .... By analogy, the Π) <sup>2</sup>*P*(*x*, *τ*) term is found. By condition *III*, matrix Δ) = <sup>6</sup> ) *P*1(*x*) *P*<sup>0</sup> Π) <sup>1</sup>*P*(*x*, *τ*) Π0*P*(*τ*) 7 has an inverse for all *x* ∈ *X*, *t* ∈ [0, *T*], *ε* ∈ (0,*ε*0], *τ* = −1, −2, ... , 1 − *T* and it follows that *<sup>Y</sup>*)−<sup>1</sup> = *<sup>E</sup>* <sup>−</sup>*Y*12*Y*22−<sup>1</sup> <sup>0</sup> <sup>Δ</sup>)−<sup>1</sup> exists and system (19), (20) is uniquely solvable with a solution ⎛ ⎜⎜⎝ M1(*x*) ΠM1(*x*, *τ*) *N*1(*x*) *N*2(*x*) ⎞ ⎟⎟⎠ <sup>=</sup> *<sup>Y</sup>*)−<sup>1</sup> ⎛ ⎜⎜⎝ *P*1(*x*) 0 −*P*2(*x*) −Π2*P*(*x*, *τ*) ⎞ ⎟⎟⎠.

Thus, the remaining coefficients of the Padé approximation M1(*x*), ΠM1(*x*, *τ*), *N*1(*x*), *N*2(*x*) are found in the form

$$\begin{array}{ll} \left( \begin{array}{c} N\_{1} \\ N\_{2} \end{array} \right) = \breve{\Lambda}^{-1} \left( \begin{array}{c} -\breve{\mathsf{P}}\_{2}(\textbf{x}) \\ -\breve{\Pi}\_{2}\mathrm{P}(\textbf{x},\tau) \end{array} \right), \\\left( \begin{array}{c} \mathbf{M}\_{1}(\textbf{x}) \\ \Pi\mathbf{M}\_{1}(\textbf{x},\tau) \end{array} \right) = \left( \begin{array}{c} \breve{\mathsf{P}}\_{1}(\textbf{x}) \\ 0 \end{array} \right) - \left( \begin{array}{c} -\breve{\mathsf{P}}\_{0} & 0 \\ -\Pi\_{0}\mathrm{P}(\tau) & 0 \end{array} \right) \breve{\Lambda}^{-1} \left( \begin{array}{c} -\breve{\mathsf{P}}\_{2}(\textbf{x}) \\ -\breve{\Pi}\_{2}\mathrm{P}(\textbf{x},\tau) \end{array} \right). \end{array} \tag{25}$$

Under condition *IV*, there exists a corresponding one-point matrix *PA* [1/2] of form (18), and the symmetric matrix of the gain coefficients of the regulator is found from (24). 

#### **4. Computational Experiments**

One of the ways to concretize the coefficients of the system of linear equations for the Padé approximation can be associated with the analysis of the coefficients of the system, and the assumption that the Padé structures form a certain framework, within which the coefficients can be improved from the point of view of the optimality criterion of the control problem. This approach is illustrated below on an example of a simple pendulum [25].

$$\begin{aligned} & \mathbf{t} \in [0, 1], \; \varepsilon = 0.05, \; N = 20, \\ & A(\mathbf{x}) = \begin{pmatrix} 1 & 0.01 \\ \frac{-10 \sin(\mathbf{x}\_1)}{\mathbf{x}\_1} & 1 \end{pmatrix}, B(\mathbf{x}) = \begin{pmatrix} 0.1 \\ 0.1 \end{pmatrix}, A\_0 = \begin{pmatrix} 1 & 0.01 \\ -10 & 1 \end{pmatrix}, A\_1(\mathbf{x}) = 1/\varepsilon \begin{pmatrix} 0 & 0 \\ \frac{-10 \sin(\mathbf{x}\_1)}{\mathbf{x}\_1} + 10 & 0 \end{pmatrix}, \\ & \mathbf{x}\_0 = [0.1; -0.1]^T \\ & Q\_0 = \begin{pmatrix} 10 & 0 \\ 0 & 0.05 \end{pmatrix}, Q\_1 = \begin{pmatrix} 300 + \mathbf{x}\_1^2 & 0 \\ 0 & 300 + \mathbf{x}\_2^2 \end{pmatrix}, Q\_2 = \begin{pmatrix} 0 & 0 \\ 0 & 0 \end{pmatrix}, R = 1, \; F = \begin{bmatrix} 300 & 0 \\ 0 & 300 \end{bmatrix}. \end{aligned}$$

Here, *x*<sup>1</sup> denotes the pendulum angle, and *x*<sup>2</sup> is the angular velocity.

Using this example, taking into account the fulfillment of the conditions for the existence of matrices in the representations (18), we will demonstrate the Algorithm 1 for discrete modified Padé regulator construction and the results of its work.

In Table 1, the criterion values for the regulators based on the uniform first-order asymptotic approximation (*P*0, ) *P*1(*x*), M0*P*(*τ*), Π) <sup>1</sup>*P*(*x*, *τ*)) and the modified Padé [1/2] from (22) are presented. The comparison is made with the D-SDRE regulator which uses the solution of Equation (9) by the algorithm from [5].

**Table 1.** Regulators comparison by criterion values for *ε* = 0.05.


**Algorithm 1:** Discrete modified Padé regulator construction.



The introduction of multipliers allows us to correct the *PA* system matrix coefficients within the obtained structure. Such a technique can be used as a basis for the correction of the resulting Padé regulator if the result by the optimality criterion is better than the corresponding results along the regulators using only the asymptotics and regulators built based on the SDRE technique, which is demonstrated in the calculations given below. A regulator built based on this approach will be called a modified Padé regulator.

4. The resultant modified Padé regulator gain is found from (23).

Table 2 shows that the modified Padé [1/2] regulator is closer to the D-SDRE solution by optimality criterion values on a larger interval of parameter variation in comparison with the regulator based on the uniform first-order asymptotic expansion, i.e., demonstrates good quality approximation for larger values of the parameter and has better extrapolation properties. In this example, the uniform first-order asymptotic approximation works only for small parameter values and fails to stabilize the system when the value of the parameter increases. Moreover, the Padé regulator is significantly better by criterion value than the D-SDRE regulator. Thus, here the modified Padé regulator is sufficiently better by optimality criterion than the two other control algorithms (D-SDRE, asymptotic approximation) in the selected parameter variation interval from 0.05 to 0.25 and the asymptotic approximation has a restricted area of application and provides worse quality of the approximation.


**Table 2.** Demonstration of extrapolation properties of the modified Padé regulator.

The corresponding closed-loop trajectories are presented in Figure 1. It can be seen that the Modified Padé [1/2] approximation brings the system to the neighborhood of the zero-equilibrium point and the trajectories are similar to the D-SDRE solution.

**Figure 1.** Closed-loop system trajectories: Modified Padé [1/2] approximation (red lines) and D-SDRE (blue lines).

#### **5. Conclusions**

Using the SDRE approach, the asymptotics of the solution of the corresponding initial singularly perturbed control problem for the matrix discrete Riccati equation with coefficients weakly dependent on the state is constructed and the corresponding one-point *PA* regulator is proposed, i.e., only one asymptotic approximation of the Riccati equation solution is used to construct the *PA* for the feedback gain matrix of the regulator. The results of numerical experiments illustrate, in particular, the improved extrapolation properties of the constructed regulator, which makes the algorithm applicable in control systems for a wider range of parameter variation. An approach for modified *PA* construction is also demonstrated, which consists of the correction of the system of equations for finding the *PA* coefficients taking into account the structure of the matrix of the original system and the properties of the terms of the asymptotic approximation.

**Author Contributions:** Conceptualization, M.D.; methodology, M.D. and Y.D.; software, Y.D.; formal analysis, M.D.; writing—original draft preparation, M.D. and Y.D.; writing—review and editing, M.D. and Y.D. All authors have read and agreed to the published version of the manuscript.

**Funding:** Research is supported by the Russian Science Foundation (Project No. 21-11-00202).

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Not applicable.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**

