Optimization of Weight Matrices for the Linear Quadratic Regulator Problem Using Algebraic Closed-Form Solutions

Abstract

This work proposes an analytical gradient-based optimization approach to determine the optimal weight matrices that make the state and control input at the final time close to zero for the linear quadratic regulator problem. Most existing methodologies focused on regulating the diagonal elements using only bio-inspired approaches or analytical approaches. The method proposed, contrarily, optimizes both diagonal and off-diagonal matrix elements based on the gradient. Moreover, by introducing a new variable composed of the steady-state and time-varying terms for the Riccati matrix and using the coordinate transformation for the state, one develops algebraic equationsbased closed-form solutions to generate the required states and numerical partial derivatives for an optimization strategy that does not require the computationally intensive numerical integration process. The authors test the algorithm with one- and two-degrees-of-freedom linear plant models, and it yields the weight matrices that successfully satisfy the pre-defined requirement, which is the norm of the augmented states less than 10⁻⁵. The results suggest the broad applicability of the proposed algorithm in science and engineering.

1. Introduction

In optimal control theory, one defines optimal control problems to determine control signals that minimize a performance index as well as satisfy the physical constraints [1]. Whereas it is difficult to find explicit expressions for the optimal control of nonlinear systems, linear systems can provide explicit equations for the optimal control input. For the advantage of using explicit equations, many works linearize nonlinear systems at the equilibrium point. One of the widely used state feedback control problems is defined for a linear system, where a quadratic performance index is defined using the state, control, terminal state, and associated unspecified weight matrices, called the linear quadratic regulator (LQR) problem. In general, the weight matrices are adjusted by users to find the optimal control and corresponding states that satisfy the requirement or desired behavior set by users. For instance, if it is important that the intermediate state is to be small, a large state weight matrix can be selected. The other weight matrices, such as the terminal state and control weight matrices concerning the terminal state and control effort, are selected in a similar manner. Although a key of solving the LQR problem is to determine the weight matrices properly, most studies have determined these matrices by trial-and-error, which is time-consuming work. For instance, considering a second-order system with one degree-of-freedom (DOF) requires determining seven symmetric components of the weight matrices (i.e., three for the state weight $Q\in {\mathbb{R}}^{2\times 2}$, one for the control weight $R\in \mathbb{R}$, and three for the terminal state weight $S_{\mathrm{f}}\in {\mathbb{R}}^{2\times 2}$). However, if one considers a two-DOFs system, the number of the symmetric components in the weight matrices increases to 23 (i.e., 10 for $Q\in {\mathbb{R}}^{4\times 4}$, three for $R\in {\mathbb{R}}^{2\times 2}$, and 10 for $S_{\mathrm{f}}\in {\mathbb{R}}^{4\times 4}$). That is, the number of components to be determined in the weight matrices exponentially increases if the size of the system increases. Moreover, determining the weight matrices generally depends on the engineers’ knowledge, and this approach does not always guarantee that the system response will meet the specified performance. For this reason, it is critical to design a proper optimization strategy for determining the weight matrices while satisfying the user-specified performance and reducing the engineers’ effort to find the proper weight matrices, especially for high DOF applications.

Many researchers have proposed weight matrices selection approaches to reduce the effort for determining weight matrices using bio-inspired heuristic algorithms. Among a number of bio-inspired optimization algorithms, some researchers employed a genetic algorithm (GA) to optimize the weight matrices for various system models and fitness functions. Marada et al. [2] formulated the LQR problem based on the linearized inverted pendulum system and found the state and control weight matrices using the GA. Their approach minimizes a fitness function that consists of a weighted summation of the settling time and the control input considering the location of system closed-loop poles for system stability. Dhiman et al. [3] also utilized the same system model and optimization algorithm, but the considered fitness function was the quadratic performance index. The GA-based weight optimization approach was also utilized for the vehicle suspension system by Yu et al. [4] and linearized spacecraft attitude dynamics by Kukreti et al. [5]. For the vehicle suspension system, the fitness function was composed of suspension performance indices, such as vertical acceleration of the car, dynamic tire load, and suspension working space. Similarly, Kukreti et al. considered the control performance for the fitness function, which is defined as the weighted summation of the final time and attitude and angular velocity errors. In addition to the GA, to optimize the weight matrices of the LQR problem, scholars applied many different optimization algorithms (e.g., a particle swarm optimization [6], its variations [7,8], ant-lion algorithm [9], ant colony optimization [10], Jaya’s algorithm [11], and bat algorithm [12], etc.) into diverse dynamic systems, and more studies including the mentioned research are summarized in Table 1. The aforementioned approaches commonly considered the diagonal components of the state and control weight matrices, as well as the algebraic Riccati matrix equation. Also, these approaches, based on bio-inspired optimization algorithms, take offline optimization which requires a huge computational burden.

Unlike weight selection using bio-inspired heuristic algorithms, some scholars proposed analytical approaches to determine the weight matrices. Elumalai and Raaja [13] proposed a time-domain-based algebraic Riccati and Lagrange optimization technique to determine the weight matrices, where they satisfy a pre-defined performance for a two-dimensional torsion system. Sarkar and Dewan [14] presented a pole-placement approach that satisfies user-defined time domain specifications for generating the state-feedback gain for an inverted pendulum system. Yang [15] presented a pole assignment design of a quaternion-based spacecraft attitude control problem for generating the weight matrices by considering a balance between the performance and fuel consumption.

Table 1. Relevant research papers.

Refs.	System Model	Approach	Consideration
[2]	Inverted pendulum	Genetic algorithm	- Weighted summation of settling time and control input
[3]	Inverted pendulum	Genetic algorithm	- Performance index of the LQR problem
[4]	Active suspension system	Genetic algorithm	- Vertical acceleration of the car, dynamic tire load, and suspension working space
[5]	Spacecraft attitude dynamics	Genetic algorithm	- Weighted summation of the final time, attitude, and angular velocity error
[6]	Inverted double pendulum	Particle swarm optimization	- Performance index of the LQR problem
[7]	Inverted pendulum	Adaptive particle swarm optimization	- Integrated state error
[8]	Semi-active suspension system	Improved particle swarm optimization	- Vertical acceleration, suspension deflection and tire dynamic load of the passive suspension vehicle
[9]	Inverted double pendulum	Ant-lion algorithm	- Crane position, upswing and downswing angle
[10]	Vehicle suspension system	Ant colony optimization	- Suspension travel, suspension velocity, tire deflection, and tire velocity
[11]	Deregulated power system	Jaya’s algorithm	- State error
[12]	Vehicle suspension system	Bat algorithm	- Ride comfort and passenger safety
[16]	Vehicle suspension system	Adaptive predator-prey optimization	- Integral square error of the state
[17]	Two-axis CNC system	Artificial bee colony optimization	- Settling time and overshoot
[18]	Inverted double pendulum	Neuro Evolution of Augmenting Topologies	- State error
[13]	Two-dimensional torsion system	Lagrange optimization	- System response (overshoot, setting time, steady state error)
[14]	Inverted pendulum	Pole-placement approach	- Time-domain specifications
[15]	spacecraft attitude dynamics	Pole assignment design	- Control performance and fuel consumption

Table 1. Relevant research papers.

Refs.	System Model	Approach	Consideration
[2]	Inverted pendulum	Genetic algorithm	- Weighted summation of settling time and control input
[3]	Inverted pendulum	Genetic algorithm	- Performance index of the LQR problem
[4]	Active suspension system	Genetic algorithm	- Vertical acceleration of the car, dynamic tire load, and suspension working space
[5]	Spacecraft attitude dynamics	Genetic algorithm	- Weighted summation of the final time, attitude, and angular velocity error
[6]	Inverted double pendulum	Particle swarm optimization	- Performance index of the LQR problem
[7]	Inverted pendulum	Adaptive particle swarm optimization	- Integrated state error
[8]	Semi-active suspension system	Improved particle swarm optimization	- Vertical acceleration, suspension deflection and tire dynamic load of the passive suspension vehicle
[9]	Inverted double pendulum	Ant-lion algorithm	- Crane position, upswing and downswing angle
[10]	Vehicle suspension system	Ant colony optimization	- Suspension travel, suspension velocity, tire deflection, and tire velocity
[11]	Deregulated power system	Jaya’s algorithm	- State error
[12]	Vehicle suspension system	Bat algorithm	- Ride comfort and passenger safety
[16]	Vehicle suspension system	Adaptive predator-prey optimization	- Integral square error of the state
[17]	Two-axis CNC system	Artificial bee colony optimization	- Settling time and overshoot
[18]	Inverted double pendulum	Neuro Evolution of Augmenting Topologies	- State error
[13]	Two-dimensional torsion system	Lagrange optimization	- System response (overshoot, setting time, steady state error)
[14]	Inverted pendulum	Pole-placement approach	- Time-domain specifications
[15]	spacecraft attitude dynamics	Pole assignment design	- Control performance and fuel consumption

Most of the existing studies only find diagonal components for the state and control weight matrices using bio-inspired heuristic optimization algorithms along with the infinite-time algebraic Riccati matrix equation, and the limited number of studies deals with analytical methods for finding diagonal components of the state and control weight matrices as shown in Table 1. Considering only diagonal components for the weight matrices may not provide better control performance if the states of the considered system are coupled with each other [19]. On the other hand, this work optimizes all symmetric components of the weight matrices in the LQR problem to minimize the state and control values at the final time by utilizing analytic gradients. Also, the optimization process only contains algebraic equations for the principal equations and the partial derivatives, which are developed by employing the steady-state and time-varying terms for the time-varying differential Riccati matrix equation [20] and utilizing the coordinate transformation for the states [21]. The main contributions of this work include the following:

• We design a gradient-based optimization strategy for determining diagonal and off-diagonal components of the state, control, and terminal state weight matrices; to provide more flexibility for optimization.
• We find only algebraic equations-contained closed-form solutions for the principal equations and the sensitivity partial derivatives, including the time-varying Riccati matrix equation, which require less computational cost for the optimization.

As a result, no numerical integration is required for solving the finite-time fixed-time optimal feedback control as well as optimizing the weight matrices. In simulation studies, two most widely used numerical examples are presented to demonstrate the effectiveness of the proposed approach, for second-order linear differential systems with one and two DOFs.

2. Formulation of Linear Quadratic Regulator Problem

The optimal state feedback control problem for a linear time-invariant and continuous-time model with a finite and fixed-time, called the LQR problem, is to find the control inputs that minimize the following performance index [22]:

$$\mathcal{L}=\frac{1}{2}{\mathbf{x}}^{\mathrm{T}}\left(t_{\mathrm{f}}\right)S_{\mathrm{f}}\mathbf{x}\left(t_{\mathrm{f}}\right)+\frac{1}{2}{\int }_{t_0}^{t_{\mathrm{f}}}\left[{\mathbf{x}}^{\mathrm{T}}\left(t\right)Q\mathbf{x}\left(t\right)+{\mathbf{u}}^{\mathrm{T}}\left(t\right)R\mathbf{u}\left(t\right)\right]\mathrm{d}t,$$

(1)

subject to

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

(2)

with the initial condition $\mathbf{x}\left(t_0\right)={\mathbf{x}}_0$ and terminal condition $\mathbf{x}\left(t_{\mathrm{f}}\right)={\mathbf{x}}_{\mathrm{f}}$. The optimal control is obtained by [22]:

$${\mathbf{u}}^{*}\left(t\right)=-R^{-1}{B}^{\mathrm{T}}S\left(t\right)\mathbf{x}\left(t\right),$$

(3)

where the time-varying Riccati matrix is found by integrating the following backward in time:

$$\dot{S}\left(t\right)=-S\left(t\right)A-A^{\mathrm{T}}S\left(t\right)+S\left(t\right)B{R}^{-1}{B}^{\mathrm{T}}S\left(t\right)-Q.$$

(4)

This work assumes that $(A,B)$ is controllable and $(A,D)$ is observable. Note that $D\in {\mathbb{R}}^{g\times n}$ with $g\le n$ is defined as $Q=D^{\mathrm{T}}D$ and can be obtained by Cholesky decomposition. In addition, the proof of the stability of the given closed-loop control system is discussed in Appendix A. To find the optimal control ${\mathbf{u}}^{*}\left(t\right)$, indeed, it is required to properly select the weight matrices (Q, R, and $S_{\mathrm{f}}$), which are the user-defined parameters in general. Then, one can compute the time-varying Riccati matrix, state, and control input through numerical integration. If the user-defined requirement, such as the norm of the states less than a certain value at the final time, is satisfied, one can obtain the state with the optimal control input as shown in Figure 1. In general, the weight matrices are determined by the engineers’ knowledge and trial-and-error. The engineers need to change the weight matrices iteratively until the requirement is satisfied. However, it is not only time-consuming work but also hard to guarantee that the system response satisfies the specified performance.

3. Optimization of Weight Matrices

To efficiently determine the weight matrices of the LQR problem, this work proposes an optimization process for the weight matrices. Unlike the aforementioned studies in the introduction, this work considers all symmetric components of the weight matrices (Q, R, and $S_{\mathrm{f}}$) in the formulation to gain more flexibility in optimization. Here, the unknown parameters for the optimization process are the symmetric elements of Q, R, and $S_{\mathrm{f}}$ are as follows:

$$\begin{array}{cc}\hfill Q& =\left[\begin{array}{ccc}{q}_1& \cdots & q_n\\ ⋮& \ddots & ⋮\\ q_n& \cdots & q_N\end{array}\right],\hfill \end{array}$$

(5)

$$\begin{array}{cc}\hfill R& =\left[\begin{array}{ccc}{r}_1& \cdots & r_m\\ ⋮& \ddots & ⋮\\ r_m& \cdots & r_M\end{array}\right],\hfill \end{array}$$

(6)

$$\begin{array}{cc}\hfill S_{\mathrm{f}}& =\left[\begin{array}{ccc}{s}_1& \cdots & s_n\\ ⋮& \ddots & ⋮\\ s_n& \cdots & s_N\end{array}\right].\hfill \end{array}$$

(7)

For the simple notation, the symmetric components of the weight matrices are gathered into one vector as

$$\begin{array}{cc}\hfill \mathbf{w}& \equiv {[w_1\cdots w_P]}^{\mathrm{T}}\in {\mathbb{R}}^P\hfill \\ \hfill & ={\left[q_1\cdots q_N{r}_1\cdots r_M{s}_1\cdots s_N\right]}^{\mathrm{T}}.\hfill \end{array}$$

(8)

For optimization, this work introduces a new variable consisting of the state and control, called an augmented state, as follows:

$$\mathbf{y}\left(t\right)=\left[\begin{array}{c}\mathbf{x}\left(t\right)\\ \mathbf{u}\left(t\right)\end{array}\right]\in {\mathbb{R}}^{n+m}.$$

(9)

The goal of the proposed optimization process is to find the optimal weight matrices that minimize the performance index while $\mathbf{y}\left(t\right)$ at the final time is to be zero. Thus, for Equation (9), one applies Taylor expansion at $t=t_{\mathrm{f}}$ as a function of the weight matrices’ symmetric elements defined in Equations (5)–(7), which leads to

$$\mathbf{y}\left(t_{\mathrm{f}}\right)=\mathbf{y}(t_{\mathrm{f}},Q,R,S_{\mathrm{f}})+\sum _{k=1}^N\mathbf{y}{,}_{q_k}\mathrm{d}{q}_k+\sum _{l=1}^M\mathbf{y}{,}_{r_l}\mathrm{d}{r}_l+\sum _{k=1}^N\mathbf{y}{,}_{s_k}\mathrm{d}{s}_k.$$

(10)

It is important to note that $\beta {,}_{\alpha }$ indicates the partial derivative of an arbitrary variable $\beta$ with respect to an arbitrary variable $\alpha$. That is, “,” between the variable and the subscript indicates the partial derivative. Targeting $\mathbf{y}\left(t_{\mathrm{f}}\right)=\mathbf{0}$ and collecting the partial derivatives in Equation (10) into a global Jacobian matrix J leads to

$$\begin{array}{cc}\hfill \mathbf{0}& =\mathbf{y}(t_{\mathrm{f}},Q,R,S_{\mathrm{f}})+J\mathrm{d}\mathbf{w}\hfill \\ \hfill & ={\mathbf{y}}_{\mathrm{f}}+J\mathrm{d}\mathbf{w},\hfill \end{array}$$

(11)

where

$$\begin{array}{cc}\hfill J=& \left[\mathbf{y}{,}_{q_1}\cdots \mathbf{y}{,}_{q_N}\mathbf{y}{,}_{r_1}\cdots \mathbf{y}{,}_{r_M}\mathbf{y}{,}_{s_1}\cdots \mathbf{y}{,}_{s_N}\right]\in {\mathbb{R}}^{(n+m)\times P},\hfill \end{array}$$

(12)

$$\begin{array}{cc}\hfill \mathrm{d}\mathbf{w}=& {\left[\mathrm{d}{q}_1\cdots \mathrm{d}{q}_N\mathrm{d}{r}_1\cdots \mathrm{d}{r}_M\mathrm{d}{s}_1\cdots \mathrm{d}{s}_N\right]}^{\mathrm{T}}\in {\mathbb{R}}^P.\hfill \end{array}$$

(13)

Here, each element in the correction vector $\mathrm{d}\mathbf{w}$ represents changes in each symmetric element of the weight matrices. Since the number of unknown parameters ($P=2N+M$) is larger than the number of final conditions ($n+m$), one can obtain the solution for $\mathrm{d}\mathbf{w}$ by minimizing the following:

$$\mathcal{H}=\frac{\mathrm{d}{\mathbf{w}}^{\mathrm{T}}\mathrm{d}\mathbf{w}}{2}+{\mathbf{\lambda }}^{\mathrm{T}}\left({\mathbf{y}}_{\mathrm{f}}+J\mathrm{d}\mathbf{w}\right),$$

(14)

The necessary conditions for the optimization are derived as

$$\begin{array}{cc}\hfill \mathcal{H}{,}_{\mathrm{d}\mathbf{w}}& =\mathrm{d}\mathbf{w}+J^{\mathrm{T}}\mathbf{\lambda }=\mathbf{0},\hfill \end{array}$$

(15)

$$\begin{array}{cc}\hfill \mathcal{H}{,}_{\mathbf{\lambda }}& ={\mathbf{y}}_{\mathrm{f}}+J\mathrm{d}\mathbf{w}=\mathbf{0}.\hfill \end{array}$$

(16)

Then, manipulating two necessary conditions, the solution for $\mathbf{\lambda }$ is given by

$$\mathbf{\lambda }={\left(J{J}^{\mathrm{T}}\right)}^{-1}{\mathbf{y}}_{\mathrm{f}},$$

(17)

and the minimum norm optimization solution for the parameter correction vector is obtained as

$$\mathrm{d}\mathbf{w}=-J^{\mathrm{T}}{\left(J{J}^{\mathrm{T}}\right)}^{-1}{\mathbf{y}}_{\mathrm{f}}.$$

(18)

Therefore, the symmetric components of the weight matrices are updated as

$${\mathbf{w}}_{\mathrm{update}}={\mathbf{w}}_{\mathrm{previous}}+\mathrm{d}\mathbf{w}.$$

(19)

Figure 2 displays the procedure of the proposed optimization process.

Once initial weight matrices are assumed, one needs to confirm whether $S_{\mathrm{f}}$ is equal to $S_{\mathrm{ss}}$ or not. This condition comes from deriving the closed-form solution for the time-varying Riccati matrix that will be explained in Section 4.1.1. If $S_{\mathrm{f}}=S_{\mathrm{ss}}$ at the first step, users need to redefine the initial $S_{\mathrm{f}}$. However, if $S_{\mathrm{f}}=S_{\mathrm{ss}}$ in a second or higher iteration, $S_{\mathrm{f}}$ obtained is slightly changed by replacing it with $(1-\gamma )S_{\mathrm{f}}$. Note that $\gamma$ is a user-defined small constant value, which is considered not only to avoid the equality condition but also to use the values as close to the updated values as possible. After that, one requires to verify whether the weight matrices satisfy the definiteness condition or not. That is, Q and $S_{\mathrm{f}}$ must be positive semi-definite, and R must be positive definite. In particular, after finding D matrix via Cholesky decomposition, the observability for $(A,D)$ pair is evaluated. After confirming the definiteness condition for each updated weight matrix, one finds $S\left(t\right)$, $\mathbf{x}\left(t\right)$, and $\mathbf{u}\left(t\right)$, and evaluates ${\mathbf{y}}_{\mathrm{f}}$ using the confirmed weight matrices only. Note that it uses previous weight matrices if the updated weight matrices violate the definiteness conditions. For instance, if all weight matrices violate the definiteness condition, one should try different initial weight matrices. However, for example, if only the updated Q violates the definiteness condition, it uses Q obtained from the previous step and the updated R and $S_{\mathrm{f}}$ for the next step. Next, it evaluates the two-stage stopping conditions. First, it terminates the process when the current iteration number exceeds the maximum iteration number, which is defined by users, and then repeats the entire process with the newly assumed initial weight matrices. This means that the initial guesses for the weight matrices used are not properly selected, and it requires starting from new initial guesses to find the optimal weight matrices. Otherwise, it evaluates the second stopping condition. If the resulting solution meets the requirement defined by users ($ϵ$), it terminates the process and endorses the weight matrices as optimal. If not, it updates the weight matrices. To proceed with the update process, one evaluates the Jacobian J. After computing $\mathrm{d}\mathbf{w}$ and updating $\mathbf{w}$ using $\mathrm{d}\mathbf{w}$, one converts ${\mathbf{w}}_{\mathrm{update}}$ into the corresponding updated weight matrix. This procedure continues until the resulting solution satisfies the requirement.

The proposed optimization process requires finding the state and control using Equations (2)–(4) to evaluate the augmented state at the final time and computing the partial derivatives of Equations (2)–(4) with respect to each symmetric component of the weight matrices to determine the Jacobian matrix in Equation (12). Then, the partial derivatives for the state, control, and time-varying Riccati matrix with respect to each symmetric element of the weight matrices $w_p$ (for $p=1,\dots ,P$) are derived as

$$\begin{array}{cc}\hfill \dot{\mathbf{x}}{,}_{w_p}\left(t\right)=& A\mathbf{x}{,}_{w_p}\left(t\right)+B\mathbf{u}{,}_{w_p}\left(t\right),\hfill \end{array}$$

(20)

$$\begin{array}{cc}\hfill \mathbf{u}{,}_{w_p}\left(t\right)=& R^{-1}R{,}_{w_p}{R}^{-1}{B}^{\mathrm{T}}S\left(t\right)\mathbf{x}\left(t\right)-R^{-1}{B}^{\mathrm{T}}S{,}_{w_p}\left(t\right)\mathbf{x}\left(t\right)-R^{-1}{B}^{\mathrm{T}}S\left(t\right)\mathbf{x}{,}_{w_p}\left(t\right),\hfill \end{array}$$

(21)

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill \dot{S}{,}_{w_p}\left(t\right)=& -S{,}_{w_p}\left(t\right)A-A^{\mathrm{T}}S{,}_{w_p}\left(t\right)+S{,}_{w_p}\left(t\right)B{R}^{-1}{B}^{\mathrm{T}}S\left(t\right)\hfill \\ \hfill & -S\left(t\right)B{R}^{-1}R{,}_{w_p}{R}^{-1}{B}^{\mathrm{T}}S\left(t\right)+S\left(t\right)B{R}^{-1}{B}^{\mathrm{T}}S{,}_{w_p}\left(t\right)-Q{,}_{w_p},\hfill \end{array}\end{array}$$

(22)

where $\mathbf{x}{,}_{w_p}\left(t_0\right)=\mathbf{0}$ and $S{,}_{w_p}\left(t_0\right)=0$, because $\mathbf{x}\left(t_0\right)$ and $S\left(t_0\right)$ are constants. So, to evaluate the augmented state, it first numerically integrates the differential Riccati matrix equation in Equation (4) backward in time and then integrates the controlled system response computed by substituting the optimal control defined in Equation (3) into Equation (2) forward in time. In the same manner, to compute the Jacobian matrix, it integrates the partial derivatives in Equations (20) and (22) forward in time. In fact, this process is technically straightforward. However, it requires the exponential computational load as the dimension of the state increases because of multiple numerical integrations. Hence, this work develops closed-form algebraic expressions for the principal equations and partial derivatives. The use of closed-form solutions not only eliminates the need for introducing numerical integration methods but also allows computing analytic sensitivity partial derivatives. Otherwise, it would require extensive numerical integration calculations.

4. Closed-Form Solutions for Principal Equations and Partial Derivatives

This section introduces the closed-form solutions for the principal equations and the sensitivity partial derivatives. The closed-form solutions for the principal equations are mainly used to solve the LQR problem and find the augmented state at the final time, and the closed-form solutions for the partial derivatives are utilized to obtain the Jacobian in the optimization process.

4.1. Derivation of Closed-Form Solutions for Principal Equations

This work considers the time-varying Riccati matrix, the controlled state, control input, and the performance index as the principal equations. It is important to note that the closed-form solutions derived here do not require numerical integration, so there is a minimal computational load to solve the LQR problem.

4.1.1. Time-Varying Riccati Matrix

It eliminates numerical integration by introducing the following closed-form solution for the differential Riccati matrix equation in Equation (4), where the solution consists of a steady-state term and a time-varying term [20,23]:

$$S\left(t\right)=S_{\mathrm{ss}}+Z^{-1}\left(t\right).$$

(23)

Note that $Z\left(t\right)$ is invertible, and the detailed condition is described after Equation (27). Also, the final condition is given by $S_{\mathrm{f}}$, and the steady-state solution $S_{\mathrm{ss}}$ satisfies the algebraic Riccati matrix equation:

$$-S_{\mathrm{ss}}A-A^{\mathrm{T}}{S}_{\mathrm{ss}}+S_{\mathrm{ss}}B{R}^{-1}{B}^{\mathrm{T}}{S}_{\mathrm{ss}}-Q=0.$$

(24)

Substituting Equation (23) into Equation (4), the differential Lyapunov matrix equation for $Z\left(t\right)$ is derived as

$$\dot{Z}\left(t\right)=\overline{A}Z\left(t\right)+Z\left(t\right){\overline{A}}^{\mathrm{T}}-B{R}^{-1}{B}^{\mathrm{T}},$$

(25)

where $\overline{A}$ is defined as [24]

$$\overline{A}=A-B{R}^{-1}{B}^{\mathrm{T}}{S}_{\mathrm{ss}},$$

(26)

and from Equation (23), the final condition of $Z\left(t\right)$ is found to be

$$Z\left(t_{\mathrm{f}}\right)={(S_{\mathrm{f}}-S_{\mathrm{ss}})}^{-1}.$$

(27)

Note that the condition $S_{\mathrm{f}}\ne S_{\mathrm{ss}}$ must be satisfied during the optimization process so that $S_{\mathrm{f}}-S_{\mathrm{ss}}$ is invertible. This condition aligns with the reversibility of $Z\left(t\right)$ in Equation (23). The closed-form solution for $Z\left(t\right)$ with a steady-state term and a time-varying term is derived as

$$Z\left(t\right)=Z_{\mathrm{ss}}+e^{\overline{A}(t-t_{\mathrm{f}})}{Z}_{\mathrm{b}}{e}^{\overline{A}(t-t_{\mathrm{f}})},$$

(28)

where $e^{(·)}$ is the ${\mathbb{R}}^{n\times n}$ exponential matrix, and the boundary condition of $Z\left(t\right)$ is defined by

$$Z_{\mathrm{b}}={(S_{\mathrm{f}}-S_{\mathrm{ss}})}^{-1}-Z_{\mathrm{ss}}.$$

(29)

Here, the steady-state solution for $Z_{\mathrm{ss}}$ satisfies the algebraic Lyapunov matrix equation [25]:

$$\overline{A}{Z}_{\mathrm{ss}}+Z_{\mathrm{ss}}{\overline{A}}^{\mathrm{T}}-B{R}^{-1}{B}^{\mathrm{T}}=0.$$

(30)

As described above, several equations and solutions for the algebraic Riccati matrix equation and algebraic Lyapunov matrix equation are interconnected with each other to obtain the closed-form solution for the time-varying Riccati matrix. To sum up, the time-varying Riccati matrix can be obtained by substituting the solution of Equations (24) and (28) into Equation (23) as the closed-form solution. See Appendix B.1 for derivations of Equations (25) and (28).

4.1.2. State for the Closed-Loop Control System

From the problem formulation, the governing differential equation for the controlled state is derived by substituting Equations (3) and (23) into Equation (2):

$$\dot{\mathbf{x}}\left(t\right)=\left[\overline{A}-B{R}^{-1}{B}^{\mathrm{T}}{Z}^{-1}\left(t\right)\right]\mathbf{x}\left(t\right),$$

(31)

with the initial condition ${\mathbf{x}}_0$. The closed-form solution for Equation (31) is expressed as [21,26]

$$\mathbf{x}\left(t\right)=\Phi (t,t_0){\mathbf{x}}_0.$$

(32)

Note that the state transition matrix $\Phi (t,t_0)$ has the explicit form defined as

$$\Phi (t,t_0)=Z\left(t\right)e^{-{\overline{A}}^{\mathrm{T}}(t-t_0)}{Z}^{-1}\left(t_0\right),$$

(33)

where $Z\left(t\right)$ is defined by Equation (28), and $\Phi (t,t_0)$ satisfies the following properties:

$$\begin{array}{cc}\hfill \Phi (t_2,t_0)& =\Phi (t_2,t_1)\Phi (t_1,t_0),\hfill \end{array}$$

(34)

$$\begin{array}{cc}\hfill \Phi (t_0,t_1)& ={\Phi }^{-1}(t_1,t_0),\hfill \end{array}$$

(35)

$$\begin{array}{cc}\hfill \Phi (t_0,t_0)& =I.\hfill \end{array}$$

(36)

Therefore, the closed-form solution for the state is expressed as

$$\mathbf{x}\left(t\right)=\left(Z_{\mathrm{ss}}+e^{\overline{A}(t-t_{\mathrm{f}})}{Z}_{\mathrm{b}}{e}^{\overline{A}(t-t_{\mathrm{f}})}\right)e^{-{\overline{A}}^{\mathrm{T}}(t-t_0)}{Z}^{-1}\left(t_0\right){\mathbf{x}}_0.$$

(37)

See Appendix B.2 for derivations of Equations (32) and (33).

4.1.3. Control for the Closed-Loop Control System

The control input is defined by introducing Equation (23) into Equation (3):

$${\mathbf{u}}^{*}\left(t\right)=-R^{-1}{B}^{\mathrm{T}}(S_{\mathrm{ss}}+Z^{-1}\left(t\right))\mathbf{x}\left(t\right).$$

(38)

Using Equation (32) and expanding Equation (38) lead to

$${\mathbf{u}}^{*}\left(t\right)=-R^{-1}{B}^{\mathrm{T}}{S}_{\mathrm{ss}}\Phi (t,t_0){\mathbf{x}}_0-R^{-1}{B}^{\mathrm{T}}{Z}^{-1}\left(t\right)\Phi (t,t_0){\mathbf{x}}_0,$$

(39)

and substituting Equation (33) into Equation (39) yields:

$${\mathbf{u}}^{*}\left(t\right)=-R^{-1}{B}^{\mathrm{T}}{S}_{\mathrm{ss}}Z\left(t\right)e^{-{\overline{A}}^{\mathrm{T}}(t-t_0)}{Z}^{-1}\left(t_0\right){\mathbf{x}}_0-R^{-1}{B}^{\mathrm{T}}{e}^{-{\overline{A}}^{\mathrm{T}}(t-t_0)}{Z}^{-1}\left(t_0\right){\mathbf{x}}_0.$$

(40)

Finally, the closed-form solution for the control input is expressed as

$${\mathbf{u}}^{*}\left(t\right)=-R^{-1}{B}^{\mathrm{T}}\left(S_{\mathrm{ss}}\left(Z_{\mathrm{ss}}+e^{\overline{A}(t-t_{\mathrm{f}})}{Z}_{\mathrm{b}}{e}^{\overline{A}(t-t_{\mathrm{f}})}\right)+I\right)e^{-{\overline{A}}^{\mathrm{T}}(t-t_0)}{Z}^{-1}\left(t_0\right){\mathbf{x}}_0.$$

(41)

4.1.4. Performance Index

The performance index on $[t,t_{\mathrm{f}}]$$\mathcal{L}=\frac{1}{2}{\mathbf{x}}^{\mathrm{T}}\left(t\right)S\left(t\right)\mathbf{x}\left(t\right)$ is already a closed-form solution, because the closed-form solution for $\mathbf{x}\left(t\right)$ and $S\left(t\right)$ are given by Equations (37) and (23), respectively.

Algorithm 1 summarizes the process for solving the LQR problem using the closed-form solutions. The process for finding the solution for the LQR problem using the closed-form solutions contains backward computation for obtaining the Riccati matrix from the given $S_{\mathrm{f}}$ and forward computation for obtaining the state and optimal control. Note that numerical integration processes, like Runge–Kutta methods, are not incorporated. In contrast, the conventional method requires numerical integration processes, such as backward integration for the Riccati matrix and forward integration for the controlled state.

Algorithm 1 Computation procedure for solving the LQR problem using the closed-form solutions
1:	Inputs:  A, B, $t_0$, $t_{\mathrm{f}}$, ${\mathbf{x}}_0$, ${\mathbf{x}}_{\mathrm{f}}$, Q, R, $S_{\mathrm{f}}$
2:	Outputs: $\mathbf{x}\left(t\right)$, ${\mathbf{u}}^{*}\left(t\right)$
3:	$S_{\mathrm{ss}}\leftarrow 0=-S_{\mathrm{ss}}A-A^{\mathrm{T}}{S}_{\mathrm{ss}}+S_{\mathrm{ss}}B{R}^{-1}{B}^{\mathrm{T}}{S}_{\mathrm{ss}}-Q$
4:	$\overline{A}\leftarrow A-B{R}^{-1}{B}^{\mathrm{T}}{S}_{\mathrm{ss}}$
5:	$Z_{\mathrm{ss}}\leftarrow 0=\overline{A}{Z}_{\mathrm{ss}}+Z_{\mathrm{ss}}{\overline{A}}^{\mathrm{T}}-B{R}^{-1}{B}^{\mathrm{T}}$
6:	$Z\left(t_{\mathrm{f}}\right)\leftarrow {(S_{\mathrm{f}}-S_{\mathrm{ss}})}^{-1}$
7:	for $t\leftarrow t_{\mathrm{f}}$ to $t_0$ do
8:	$Z\left(t\right)\leftarrow Z_{\mathrm{ss}}+e^{\overline{A}(t-t_{\mathrm{f}})}\left[{(S_{\mathrm{f}}-S_{\mathrm{ss}})}^{-1}-Z_{\mathrm{ss}}\right]e^{\overline{A}(t-t_{\mathrm{f}})}$
9:	$S\left(t\right)\leftarrow S_{\mathrm{ss}}+Z^{-1}\left(t\right)$
10:	end for
11:	$Z^{-1}\left(t_0\right)\leftarrow Z\left(t_0\right)\ne 0$ from line 8
12:	for $t\leftarrow t_0$ to $t_{\mathrm{f}}$ do
13:	$\mathbf{x}\left(t\right)\leftarrow Z\left(t\right)e^{-{\overline{A}}^{\mathrm{T}}(t-t_0)}{Z}^{-1}\left(t_0\right){\mathbf{x}}_0$
14:	${\mathbf{u}}^{*}\left(t\right)\leftarrow -R^{-1}{B}^{\mathrm{T}}S\left(t\right)\mathbf{x}\left(t\right)$
15:	end for

4.2. Derivation of Closed-Form Solutions for Sensitivity Partial Derivatives

The partial derivative of a function with respect to the independent variable measures the sensitivity of the function. In this section, the sensitivity partial derivatives for the aforementioned equations and variables with respect to the symmetric elements of the weight matrices are derived to form the Jacobian matrix defined in Equation (12) that is used for the proposed optimization process. Note that no numerical integration is needed to compute the partial derivatives because all of the expressions are purely algebraic.

4.2.1. Time-Varying Riccati Matrix Partial Derivatives

Based on the closed-form expressions, the partial derivatives of the Riccati matrix in Equation (23) with respect to the symmetric components of the weight matrices are given by

$$S{,}_{w_p}\left(t\right)=S_{\mathrm{ss}}{,}_{w_p}-Z^{-1}\left(t\right)Z{,}_{w_p}\left(t\right)Z^{-1}\left(t\right),$$

(42)

where $S_{\mathrm{ss}}{,}_{w_p}$ is found by taking the partial derivatives of the algebraic Riccati matrix equation in Equation (24):

$$\begin{array}{c}\hfill {\overline{A}}^{\mathrm{T}}{S}_{\mathrm{ss}}{,}_{w_p}+S_{\mathrm{ss}}{,}_{w_p}\overline{A}=-S_{\mathrm{ss}}B{R}^{-1}R{,}_{w_p}{R}^{-1}{B}^{\mathrm{T}}{S}_{\mathrm{ss}}-Q{,}_{w_p}.\end{array}$$

(43)

See Appendix C.1 for this derivation. Moreover, the detailed descriptions for $Q{,}_{w_p}$, $R{,}_{w_p}$, and $S_{\mathrm{f}}{,}_{w_p}$ are explained in the first paragraph in Appendix C. Then, $Z{,}_{w_p}\left(t\right)$ is found by taking the partial derivatives of Equation (28) as follows:

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill Z{,}_{w_p}\left(t\right)=& Z_{\mathrm{ss}}{,}_{w_p}+\left[e^{\overline{A}(t-t_{\mathrm{f}})}\right]{,}_{w_p}{Z}_{\mathrm{b}}{e}^{\overline{A}(t-t_{\mathrm{f}})}\hfill \\ \hfill & +e^{\overline{A}(t-t_{\mathrm{f}})}{Z}_{\mathrm{b}}{,}_{w_p}{e}^{{\overline{A}}^{\mathrm{T}}(t-t_{\mathrm{f}})}+e^{\overline{A}(t-t_{\mathrm{f}})}{Z}_{\mathrm{b}}\left[e^{{\overline{A}}^{\mathrm{T}}(t-t_{\mathrm{f}})}\right]{,}_{w_p}.\hfill \end{array}\end{array}$$

(44)

Here, $[·]{,}_{w_p}$ denotes the partial derivative of the matrix exponential, and $Z_{\mathrm{b}}{,}_{w_p}$ is given by taking the partial derivatives of Equation (29) as follows:

$$\begin{array}{cc}\hfill Z_{\mathrm{b}}{,}_{w_p}=& -{(S_{\mathrm{f}}-S_{\mathrm{ss}})}^{-1}(S_{\mathrm{f}}{,}_{w_p}-S_{\mathrm{ss}}{,}_{w_p}){(S_{\mathrm{f}}-S_{\mathrm{ss}})}^{-1}-Z_{\mathrm{ss}}{,}_{w_p}.\hfill \end{array}$$

(45)

Also, $Z_{\mathrm{ss}}{,}_{w_p}$ is found by taking the partial derivatives of the steady-state Lyapunov matrix in Equation (30) as follows:

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill \overline{A}{Z}_{\mathrm{ss}}{,}_{w_p}+Z_{\mathrm{ss}}{,}_{w_p}{\overline{A}}^{\mathrm{T}}=& -\overline{A}{,}_{w_p}{Z}_{\mathrm{ss}}-Z_{\mathrm{ss}}\overline{A}{,}_{w_p}^{\mathrm{T}}-B{R}^{-1}R{,}_{w_p}{R}^{-1}{B}^{\mathrm{T}},\hfill \end{array}\end{array}$$

(46)

where $\overline{A}{,}_{w_p}$ is found by taking the partial derivatives of the closed-loop system dynamics matrix in Equation (26) as follows:

$$\overline{A}{,}_{w_p}=B{R}^{-1}R{,}_{w_p}{R}^{-1}{B}^{\mathrm{T}}{S}_{\mathrm{ss}}-B{R}^{-1}{B}^{\mathrm{T}}{S}_{\mathrm{ss}}{,}_{w_p}.$$

(47)

The details of the derivatives for all variables and matrices with respect to the symmetric elements of each weight matrix are explained throughout Appendix C.1, Appendix C.2, Appendix C.3, Appendix C.4 and Appendix C.5.

4.2.2. State Partial Derivatives

The state partial derivatives are generated by modeling the state as Equation (32), leading to the state partial derivatives given by

$$\mathbf{x}{,}_{w_p}\left(t\right)=\Phi {,}_{w_p}(t,t_0){\mathbf{x}}_0,$$

(48)

where $\Phi {,}_{w_p}(t,t_0)$ is expressed as

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill \Phi {,}_{w_p}(t,t_0)=& Z{,}_{w_p}\left(t\right)e^{-\overline{A}(t-t_0)}{Z}^{-1}\left(t_0\right)+Z\left(t\right)\left[e^{-\overline{A}(t-t_0)}\right]{,}_{w_p}{Z}^{-1}\left(t_0\right)\hfill \\ \hfill & -Z\left(t\right)e^{-\overline{A}(t-t_0)}{Z}^{-1}\left(t_0\right)Z{,}_{w_p}\left(t_0\right)Z^{-1}\left(t_0\right).\hfill \end{array}\end{array}$$

(49)

The detailed expressions of the partial derivative of the matrix exponential are explained in Appendix C.6.

4.2.3. Control Partial Derivatives

The partial derivatives of the optimal control in Equation (40) are expressed as

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill \mathbf{u}{,}_{w_p}\left(t\right)=& R^{-1}R{,}_{w_p}{R}^{-1}{B}^{\mathrm{T}}{S}_{\mathrm{ss}}\Phi (t,t_0){\mathbf{x}}_0-R^{-1}{B}^{\mathrm{T}}{S}_{\mathrm{ss}}{,}_{w_p}\Phi (t,t_0){\mathbf{x}}_0\hfill \\ \hfill & -R^{-1}{B}^{\mathrm{T}}{S}_{\mathrm{ss}}\Phi {,}_{w_p}(t,t_0){\mathbf{x}}_0+R^{-1}R{,}_{w_p}{R}^{-1}{B}^{\mathrm{T}}{e}^{-{\overline{A}}^{\mathrm{T}}(t-t_0)}{Z}^{-1}\left(t_0\right){\mathbf{x}}_0\hfill \\ \hfill & -R^{-1}{B}^{\mathrm{T}}\left[e^{-{\overline{A}}^{\mathrm{T}}(t-t_0)}\right]{,}_{w_p}{Z}^{-1}\left(t_0\right){\mathbf{x}}_0\hfill \\ \hfill & +R^{-1}{B}^{\mathrm{T}}{e}^{-{\overline{A}}^{\mathrm{T}}(t-t_0)}{Z}^{-1}\left(t_0\right)Z{,}_{w_p}\left(t_0\right)Z^{-1}\left(t_0\right){\mathbf{x}}_0.\hfill \end{array}\end{array}$$

(50)

The detailed derivation process is described in Appendix C.7. Note that Equations (48) and (50) are primarily utilized to form the Jacobian matrix, which is composed of the augmented state.

5. Simulation Study

5.1. Problem Descriptions

This work considers two of the most widely used example problems for the second-order differential equation with one and two DOFs to validate the efficacy of the closed-form algebraic expressions and the performance of the proposed optimization process. The system of the example problems is shown in Figure 3. The formulations described in (1) and (2) are used, and the system dynamics and control influence matrices for each problem are defined in

Table 2. Model parameters [27].

	1 DOF	2 DOFs
State Variables	$\mathbf{x}\left(t\right)={\left[x_1{x}_2\right]}^{\mathrm{T}}\in {\mathbb{R}}^2$, $u\left(t\right)\in \mathbb{R}$	$\mathbf{x}\left(t\right)={\left[x_1{x}_2{x}_3{x}_4\right]}^{\mathrm{T}}\in {\mathbb{R}}^4$, $\mathbf{u}\left(t\right)={\left[u_1{u}_2\right]}^{\mathrm{T}}\in {\mathbb{R}}^2$
Models	$A=\left[\begin{array}{cc}0& 1\\ -\frac{k_1}{m_1}& -\frac{c_1}{m_1}\end{array}\right]$, $B=\left[\begin{array}{c}0\\ -\frac{1}{m_1}\end{array}\right]$	$A=\left[\begin{array}{cccc}0& 0& 1& 0\\ 0& 0& 0& 1\\ -\frac{k_2+k_3}{m_2}& -\frac{k_3}{m_2}& -\frac{c_2+c_3}{m_2}& \frac{c_3}{m_2}\\ \frac{k_3}{m_3}& -\frac{k_3}{m_3}& \frac{c_3}{m_3}& -\frac{c_3}{m_3}\end{array}\right]$, $B=\left[\begin{array}{cc}0& 0\\ 0& 0\\ \frac{1}{m_2}& 0\\ 0& \frac{1}{m_3}\end{array}\right]$

. Also, the simulation parameters for each variable are tabulated in

Table 3. Simulation parameters [27].

Parameter	Value
	1 DOF	2 DOFs
Mass (kg)	$m_1=1$	$m_2=1$
		$m_3=1$
Spring coefficient (N/m)	$k_1=0.64$	$k_2=1$
		$k_3=0.5$
Damping coefficient (Ns/m)	$c_1=0.16$	$c_2=0.1$
		$c_3=0.1$
Initial condition $\mathbf{x}\left(t_0\right)$	${\left[1010\right]}^{\mathrm{T}}$	${\left[10100\right]}^{\mathrm{T}}$
State weight matrix Q	$I_{2\times 2}$	$I_{4\times 4}$
Control weight matrix R	1	$I_{2\times 2}$
Terminal state weight matrix $S_{\mathrm{f}}$	$I_{2\times 2}$	$I_{4\times 4}$
Final time $t_{\mathrm{f}}$ (s)	10
Time interval $\mathrm{d}t$ (s)	0.01

.

5.2. Verification of the Closed-Form Solutions for the LQR Problem

To highlight the efficacy of the closed-form solutions, especially for the principal equations, numerical simulations are conducted using the closed-form solutions that eliminate the numerical integration process. The resulting state and control trajectories computed using the closed-form solutions are compared with the results obtained by a conventional approach, which is to find the Riccati matrix by numerically integrating backward in time and the state and control trajectories by integrating forward in time. Figure 4a,b depicts the comparison results for the state and control trajectories using the parameters listed in Table 2 and Table 3. Moreover, the state and optimal control trajectories are displayed in the figures for the changes in state and optimal control trajectories over iterations with the legend of “Initial” in Section 5.3. For Figure 4, it is observed that the state and control differences are less than $5\times {10}^{-5}$ for both cases. That is, it is confirmed that there is not much difference in terms of the solution trajectories between the conventional and proposed approaches. In addition, it is proved that the proposed approach is computationally more efficient than the conventional approach because, during 1000 simulations, the average computational time using the closed-form solutions is reduced by half compared to the one using the conventional approach, as listed in

Table 4. Comparison of average computational time.

	1 DOF	2 DOFs
Conventional: numerical integrations (s)	0.118	0.129
Proposed: closed-form algebraic equations (s)	0.057	0.060

. The improvement in the computational speed for solving the LQR problem can contribute to the increase in the control frequency. In fact, the reduction in the computational time is obvious because the proposed approach does not require any numerical integration process. As the dimension of the state and control increases, the computational efficiency also increases exponentially. In addition, the closed-form solution provides a more accurate result compared to the conventional one because the closed-form solution does not contain numerical errors that occurred from numerical integration process.

5.3. Weight Matrices Optimization

To validate the performance of the proposed optimization approach with the developed closed-form solutions for the partial derivatives and principal equations, the weight matrices given in Table 3 are optimized. The optimization process aims for finding the optimal weight matrices that make the state and the control input at the final time close to zero while minimizing the corrections for the weight matrices. The requirement is set to be the norm of the augmented state at the final time that is less than ${10}^{-5}$ (i.e., $ϵ={10}^{-5}$), and the maximum number of iterations is 100. With these requirements, two example problems are considered, and the weight matrices listed in Table 3 are used as the initial weight matrices.

Table 5. Optimization results for the 1-DOF problem.

	Initial	Optimized
Norm of the augmented state $\left|\right|{\mathbf{y}}_{\mathrm{f}}\left|\right|$	$2.68\times {10}^{-2}$	$3.82\times {10}^{-6}$
States at the final time $\mathbf{x}\left(t_{\mathrm{f}}\right)$	${[1.20-1.70]}^{\mathrm{T}}\times {10}^{-2}$	${[3.28-1.94]}^{\mathrm{T}}\times {10}^{-6}$
Control at the final time $\mathbf{u}\left(t_{\mathrm{f}}\right)$	$1.70\times {10}^{-2}$	$-2.71\times {10}^{-7}$
State weight matrix Q	$I_{2\times 2}$	$\left[\begin{array}{cc}1.63& 0.08\\ 0.08& 1.02\end{array}\right]$
Control weight matrix R	1	0.19
Terminal state weight matrix $S_{\mathrm{f}}$	$I_{2\times 2}$	$\left[\begin{array}{cc}1.08& 0.31\\ 0.31& 0.50\end{array}\right]$
Performance index $\mathcal{L}$	$2.01\times {10}^2$	$1.58\times {10}^2$

shows the results obtained using the initial and the optimized weight matrices for the one-DOF example problem. The norm of the augmented state at the final time using initial weight matrices does not meet the requirement, requiring a tuning process for the weight matrices. After applying the optimization process proposed with the closed-form solutions, the norm of the augmented state of $3.82\times {10}^{-6}$ that satisfies the requirement is obtained. The history of the norm of the augmented state versus the performance index for the one-DOF example is shown in Figure 5, and it is shown that both values evaluated are decreased during optimization from the initial evaluation using the initial weight matrices. Moreover, as shown in Figure 6, Figure 7 and Figure 8, the proposed approach optimizes the off-diagonal elements of the weight matrices as well as the diagonal elements although the identity matrices are used as the initial conditions. Note that the elements in the green-dotted boxes are the same as the elements on the other side because the weight matrices are symmetric. The proposed optimization approach provides more flexibility compared to existing studies because the proposed one optimizes all components of the weight matrices. In addition, since the proposed optimization process contains the step that checks the definiteness conditions of the obtained matrices at every iteration, one guarantees that all weight matrices updated are positive (semi-)definite. The state and control trajectories for all iterations are shown in Figure 9 and Figure 10. The responses of the state and control input become faster over the iterations to satisfy the requirement. Hence, the trajectories with blue color, which use the optimized weight matrices, are converged to the target states faster than the others.

The optimization results for the two-DOFs example problem are tabulated in

Table 6. Optimization results for the 2-DOFs problem.

	Initial	Optimized
Norm of the augmented state $\left|\right|{\mathbf{y}}_{\mathrm{f}}\left|\right|$	$2.04\times {10}^{-2}$	$7.71\times {10}^{-6}$
States at the final time $\mathbf{x}\left(t_{\mathrm{f}}\right)$	$\left[\begin{array}{c}1.94\times {10}^{-2}\\ -5.26\times {10}^{-3}\\ 2.48\times {10}^{-3}\\ -6.43\times {10}^{-4}\end{array}\right]$	$\left[\begin{array}{c}4.29\\ -2.36\\ 1.58\\ 1.32\end{array}\right]\times {10}^{-6}$
Control at the final time $\mathbf{u}\left(t_{\mathrm{f}}\right)$	$\left[\begin{array}{c}-2.48\\ 0.64\end{array}\right]\times {10}^{-3}$	$\left[\begin{array}{c}-4.84\\ -2.80\end{array}\right]\times {10}^{-6}$
State weight matrix Q	$I_{4\times 4}$	$\left[\begin{array}{cccc}1.54& -0.07& 0.17& -0.04\\ -0.07& 1.56& 0.14& -0.03\\ 0.17& 0.14& 0.75& 0.26\\ -0.04& -0.03& 0.26& 0.82\end{array}\right]$
Control weight matrix R	$I_{2\times 2}$	$\left[\begin{array}{cc}0.14& 0.03\\ 0.03& 0.26\end{array}\right]$
Terminal state weight matrix $S_{\mathrm{f}}$	$I_{4\times 4}$	$\left[\begin{array}{cccc}1.17& 0.10& -0.15& -0.07\\ 0.10& 0.97& 0.06& 0.02\\ -0.15& 0.06& 0.97& -0.03\\ -0.07& 0.02& -0.03& 0.99\end{array}\right]$
Performance index $\mathcal{L}$	$2.01\times {10}^2$	$1.58\times {10}^2$

and displayed in Figure 11, Figure 12, Figure 13, Figure 14, Figure 15 and Figure 16. The proposed approach successfully finds the weight matrices with $\left|\right|{\mathbf{y}}_{\mathrm{f}}\left|\right|$ of $7.71\times {10}^{-6}$ that satisfy the predefined requirement as shown in Table 6. The history of the norm of the augmented state and the performance index over iterations is depicted in Figure 11, and the result shows that both values continuously decrease until the requirement is satisfied. The optimization results for the weight matrices are shown in Table 6, and their histories during optimization are shown in Figure 12, Figure 13 and Figure 14. It is shown that all symmetric components including the off-diagonal terms are optimized to find the optimal solution. In addition, all weight matrices over iterations are surely positive (semi-)definite because the proposed approach updates the matrix that meets the definiteness conditions. As shown in Figure 15 and Figure 16, the controlled states obtained using the optimized weight matrices are converged into zero values faster compared to the states obtained using the initial weight matrices in order to satisfy the requirements.

For the two cases, even though the identity matrices are used as the initial conditions, the optimal weight matrices that satisfy the predefined requirement are obtained within 10 iterations. Over the iterations, the updated matrices do not violate the definiteness conditions while optimizing all symmetric components of the weight matrices using the developed closed-form solutions.

Moreover, additional simulations using different initial conditions (weight matrices) for better understanding are displayed, and the optimization results are discussed in Appendix D.

6. Conclusions

This study proposes a gradient-based optimization approach to determine the symmetric components of the weight matrices for the linear quadratic regulator problem by applying Taylor’s expansion to the state and control at the final time and finding the minimum norm optimization solution. To prevent the increase in the computational burden that arises from the increase in state dimensions and multiple numerical integration steps in the optimization process, this work develops the algebraic equations that exploit the closed-form solutions for the principal equations and their partial derivatives. Through numerical simulations for the one- and two-degrees-of-freedom second-order dynamic systems, it is validated that the proposed optimization approach finds all symmetric elements of the state, control, and terminal state weight matrices that satisfy the requirement (the norm of the augmented state of less than ${10}^{-5}$) without violating the definiteness condition. Moreover, it is confirmed that the use of closed-form solutions reduces the computation time by half compared to the use of numerical integrations. In the future, a linear time-varying system model will be considered to extend the applicable systems.