Cooperative control of large ﬂexible space structure by two planar robots

This study discusses a cooperative control problem of a ﬂexible beam by two planar space robots. The goal is to control the position and orientation of the ﬂexible structure and robots, and to suppress the vibration of the ﬂexible structure. To solve this problem, a boundary controller based on a dynamic model represented by a hybrid PDE-ODE model is presented. In particular, a boundary controller using the Lyapunov method is derived, and it is showed that the system becomes Lyapunov stable and that the system in the neighbourhood of the desired conﬁguration becomes asymptotically stable. As the system has an actuation redundancy, distribute control inputs among actuators considering each actuator’s performance can be distributed. Finally, several numerical simulations are carried out to investigate the validity of the proposed boundary controller.


INTRODUCTION
In recent years, with the development of the space industry, there is a demand for techniques to perform on-orbit servicing such as construction, refuelling, inspection, repair, and debris removal. Several tasks are currently done by extravehicular activities by astronauts [1]. However, due to the influence of zerogravity, radiation etc. in outer space, it is not suitable for human beings to work for a long time. Furthermore, the available work is also restricted from the upper limit of the force and torque which human beings can demonstrate. For these reasons, it is crucial to realise the on-orbit servicing by the space robot, and various research works have been done so far [2,3]. Among orbit services expected in the future, construction and debris removal [4] are particularly important. Large space structures such as space solar power plants are considered to be several kilometres in length, and assembly on-orbit is essential [5,6]. In carrying these tasks, it is essential to have the ability to transport large components, for example, when transporting structure assembled in low orbits to high orbits. A large structure cannot be transported by a single robot is preferably transported by a plurality of robots equipped with thrusters. Research on coordinated control such as formation flight with multiple satellites has been actively made in recent years [7], but much less work has been done on cooperative control considering the This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited. © 2020 The Authors. IET Control Theory & Applications published by John Wiley & Sons Ltd on behalf of The Institution of Engineering and Technology dynamic interference of multiple robots and structures. In this research, we consider a system that grasps a structure with multiple space robots and deals with a cooperative control method.
One of the problems in transporting the large space structure is the flexibility of the structure. When putting structures on the orbit, the weight reduction of structures is actively carried out due to the restrictions on the loading amount of the launch vehicle and the launch cost. As a result, it is necessary to consider the large space structure as a flexible component [8]. Because of its flexibility, elastic vibrations that are continuously distributed in the structure are generated, so the dynamics of the system in which the space robot grasps the flexible structure is constructed by partial differential equation (PDE) and ordinary differential equation (ODE), that is the model is described as a PDE-ODE hybrid infinite-dimensional model. However, the previous studies on transportation of space structure with robots do not consider the flexibility of the structure [9][10][11] or truncate to a finite-dimensional model even if it considered [12][13][14][15]. Finite-dimensional approximated model is derived by neglecting higher order vibration modes. Therefore, the following drawbacks may occur: the high order of a controller equivalent to the order of a finite-dimensional approximation of the structure makes its implementation unsuitable in the space mission, and spillover occurred by the disregarded high-frequency characteristic makes a system unstable. From this point of view, FIGURE 1 Grasped flexible beam and two planar space robots with onelink rigid arm the controller design based on the original distributed parameter systems is proposed and the strict stability is shown for many flexible structures [16][17][18][19][20][21][22][23][24][25][26][27][28]. However, to the best of our knowledge, research on cooperative transportation of space structure by applying a controller based on an infinite-dimensional model has not been reported.
In this paper, we propose a method for cooperative transportation of flexible space structure by two space robots. The flexible space structure is modelled as a flexible beam, and we use the Euler-Bernoulli beam model for the flexible beam. The system is described in the non-linear PDE-ODE hybrid model derived by Hamilton's principle. We derive a controller based on the total energy of the system, which makes the closed-loop system globally Lyapunov stable. Furthermore, by taking advantage of the actuator redundancy of the system, we weight the burden of the input of each actuator and determine the control input. This is useful in the operation of cooperative transportation. It is crucial to design the controller considering the remaining amount of the thruster fuel and battery and the performance of each robot. Next, we prove that the proposed controller asymptotically stabilises the closed-loop system at the vicinity of the desired states of the system by using LaSalle's invariance principle, which is extended to infinite-dimensional systems. The effectiveness of the proposed controller, is supported by the numerical simulation results that the control objective is achieved by the proposed controller.
As we have already mentioned, there are many studies about the flexible planar robot based on the hybrid PDE-ODE model [16][17][18][19][20][21][22][23][24][25][26][27][28]. However, many of these studies address only the vibration suppression problem. These studies are not for cooperative transportation of flexible structure as in this paper. Therefore, unlike the other control methods for the flexible planar robot, our controller can realise the followings in addition to vibration suppression. The flexible structure is converged to the desired position and orientation, and the angles of the space robots are converged to the desired angles. Further, we design a controller using the actuator redundancy, while basing on an infinite-dimensional model. This is also a different point from the previous research about the flexible planar robot. In addition, we have recently proposed a boundary controller for a onelink flexible arm [29]. In [29], we discussed a contact-force con-trol problem of a one-link flexible arm. However, our previous study [29] considers a controlled system represented by a linear PDE model instead of the non-linear PDE-ODE model, as in this study. Further, cooperative transportation control of space structure was not considered.
The paper is outlined as follows. In Section 2, the mathematical model of the system is derived, and the design of control input, including utilisation of drive redundancy, is proposed. Section 3 presents the semi-group setting of the linearised closed-loop system, and the asymptotic stability is proven. The numerical simulation results to show the validity of the proposed method are demonstrated in Section 4. Finally, Section 5 presents our conclusions.

2.1
Controlled system Figure 1 illustrates a controlled system. The system consists of two space robots with a one-link rigid arm and a grasped flexible beam moving in the XY plane in Figure 1. The flexible beam having length L, uniform linear density , and uniform flexural rigidity EI , satisfies the Euler-Bernoulli beam hypothesis. In Figure 1, O − XY is an absolute coordinate system, and o − xy is a local coordinate system, which fixed to one end of the flexible beam grasped by Robot 1 and rotate with it. Let (r X (t ), r Y (t )) be the position vector of the local coordinate system's origin o, and (t ) be the angle between local and absolute coordinate. w(x, t ) is the transverse displacement of the flexible beam, where x is the spatial point on the local coordinate system o − xy, and t is a time. For Robot i, i = 1, 2, rigid arm rigidly grasps the flexible beam, that is, the rigid arm of each robot is rigidly connected to the flexible beam. Thus, the gripping angle q ib is constant, and let q i1 (t ) be the angle between Robot i and Link i1. In addition, F ix (t ) and F iy (t ) are translation forces by thrusters on Robot i's base. i (t ) is a torque by thrusters and Control Moment Gyro (CMG) mounted on Robot i's base, which is act on Robot i's base Center of Gravity (COG), and i1 (t ) is a joint torque of rigid arm. In this paper, we assume that the thrusters on Robot i's base can generate a thrust of arbitrary size. When using the thrusters that can generate ON/OFF thrust only, PWM control to approximate the required thrusts [30,31] are necessary, or the Lyapunov thruster selection approach [32] is needed. For such a situation, we investigate the case where the thrusters can generate on/off thrust only through the numerical simulations. The definition of other parameters are shown in Table. 1. For simplicity of notations, a dot{ } and a prime { } ′ denote time and spatial derivatives, respectively. For example, for (t ), which is a function of t ,̇(t ) = d (t )∕dt . For w(x, t ), which is a function of t and x,ẇ(x, t ) = w(x, t )∕ t and w ′ (x, t ) = w(x, t )∕ x, which are the partial derivative. Here, we assume that the effect of orbital mechanics is negligible because the operation time of the robots is sufficiently smaller than the orbital period [12]. Distance between tip of the beam and link i1's COG Then, the position vector p of a point p identified by the distance x along the flexible beam relative to the absolute coordinate system is written as ] . (1) In addition, r i , r i1 , i = 1, 2, which each denote the position vector of robot i's base COG and link i1's COG relative to O − XY , can be expressed as where Hamilton's principle is applied to derive the equations of motion. The kinetic energy E k and potential energy E p of the system are expressed as follows: where the kinetic energy of flexible beam E k B and the kinetic energy of robots E k R are as follows: ) .
Further, the virtual work is given by where r i = [r iX , r iY ] T . Thus, the following equations of motion can be obtained: with the boundary conditions is inertia matrix which is a function of , and is nonlinear terms which is a function of anḋ. For 11 (t ), 2 (t ), 21 (t )] T ∈ ℝ 8 , the following relation holds: where P( ) ∈ ℝ 7×8 is always full rank. Thus, there are many combinations of actuator inputs u that can generate the required forces and torques Q which actually exerted onto the system. This is said that the system has actuator redundancy. We believe that, to the best of our knowledge, the mathematical model of the system is based on the infinite-dimensional model, and the equations of motion of this system are the first proposed.

Actuator redundancy
The system has an actuator redundancy, and thus, we can distribute control inputs among actuators. Now, we design the control input u considering each actuator performance based on the procedures in [33,34]. Let us define the control input u as where P † = P T (PP T ) −1 is a pseudo inverse of P, N ∈ ℝ 8 is a null-space vector which satisfies N ∈ Ker(P), ∈ ℝ is a design parameter. The second term in Equation (21) does not influence the system dynamics. Non-zero and bounded vector N is uniquely determined as Equation (22). For a given actuator configuration, the goal is to find the parameter such that minimises the following cost function J weighting of each input: where W 1 is a positive definite diagonal matrix which weights each input. A necessary condition for a minimum of J with respect to the parameter is The parameter satisfies this equation is then This is minimum solution of Equation (23) because the solution satisfies 2 J ∕ 2 = N T W 1 N > 0. Hence, the optimal input vector u for minimising the cost function can be expressed as By designing W 1 , it is possible to design the input so as to reduce inputs of one robot taking into account the remaining amount of fuel and the battery. Also, it is possible to weight inputs by the CMG and the joint drive motor to be preferentially used rather than inputs by thrusters.

Boundary controller
The control objective is to construct a boundary controller satisfying as t → ∞, where r Xd , r Yd , d , q 11d , and q 21d are constant desired values.
We derive a controller such that the system becomes Lyapunov stable by using the Lyapunov method. Motivated by the total system energy, we define the following candidate Lyapunov function: wherẽ= In addition is a diagonal matrix consisting of positive constant gains. We have added the last term in Equation (26) as an artificial potential energy to insure that the desired final states is a unique minimum of V . Now, we show the positive definiteness of V (q,q). For this, we introduce the following norm: where ‖ ⋅ ‖ 2 L 2 = ∫ L 0 (⋅) 2 dx. This is an induced norm on the following Hilbert space: From the Cauchy-Schwartz inequality, it is easy to see that w 2 (L, t ) ≤ 1 ‖w ′′ (x, t )‖ 2 L 2 and {w ′ (L, t )} 2 ≤ 2 ‖w ′′ (x, t )‖ 2 L 2 , where 1 and 2 are some positive constants. Therefore, the estimation E P +̃TK p̃∕ 2 ≥ 3 [‖w ′′ (x, t )‖ 2 L 2 +̃T̃] holds for a positive constant 3 . Thus, we can obtain the positive definiteness of V (q,q) if there exists a positive constant 4 such that Here, note that E k can be rewritten as follows: whereẼ k = 0 iffq = 0, and M L is a positive definite matrix. Therefore, we obtain Equation (30), and thus, V (q,q) is positive definite. By differentiating Equation (26) with respect to t and employing Equations (11)- (19), we obtain In the derivation of this equation, note the followings: d dt where D ∈ ℝ 7 is a nonlinear vector. In Equation (33), if we set the vector Q as where is a diagonal matrix consisting of positive constant gains, then substituting Equation (36) into Equation (33) yields dV dt (q,q) = −̇̃TK ḋ̃≤ 0.
Thus, we find function V is a Lyapunov function that guarantees stability of the solution of the closed-loop system (11)- (19). Actuator inputs u can be calculated by substituting Equations (36) to (24). In this study, we assume states r X (t ), r Y (t ), (t ), q 11 (t ), q 21 (t ), w(L, t ) and w ′ (L, t ) can be measured without time delay. In actual implementation, q 11 (t ), q 21 (t ) are measured by encoders mount on joints, and to measure r X (t ), r Y (t ), (t ), w(L, t ) and w ′ (L, t ), we should sense robot i's position and attitude. This can be achieved by observer robots that mount remote vision sensors [12].

STABILITY ANALYSIS AT THE NEIGHBOURHOOD OF THE DESIRED CONFIGURATION
In this section, we discuss the asymptotic stability of the closed-loop system. However, the original non-linear PDE-ODE hybrid system is complicated and challenging to prove global asymptotic stability. That is why we consider the closedloop system at the neighbourhood of the desired state and prove the asymptotic stability.

Equations of motion in Hilbert space
We re-derive the linearised equations of motion of the closedloop system assuming that the error vector between the current state and the desired state is small. Then, we obtain the following linearised PDEs: where M s ∈ ℝ 7×7 is positive and symmetric constant inertia matrix. We show the asymptotic stability of the closed-loop system using LaSalle's invariance principle, which is extended to infinite-dimensional systems [20]. To use this principle, we need to formulate the closed-loop system to a first-order evolution equation in an appropriate Hilbert space. For this, first, we introduce the new variable This variable contributes to the movement of the equilibrium point towards its origin. Next, we introduce the following statespace H : where = [ 1 , 2 , … , 5 ] T and = [ 1 , 2 , … , 7 ] T . In the state space H , we define the following inner product: where K E = diag(k X , k Y , k , k 11 , k 21 ). It can be shown that H with the inner product (43) becomes a Hilbert space. Finally, we define the unbounded linear operator A : D(A) ⊂ H → H as follows: with domain Then, the closed-loop system can be written as the following first order evolution equation on H : where z = [y,r X ,r Y ,̃,q 11 ,q 21 ,̇y,̇r X ,̇r Y ,̇̃,q 11 ,̇q 21 Now, we obtain the following lemma as the properties of the closed-loop system (46).

Lemma 1. The operator A generates a C 0 -semi-group of contractions.
Furthermore, the operator A −1 is compact.
Proof. First, we show that the operator A is dissipative. From the simple calculations using integration by parts and the boundary conditions in (45), we obtain the following: for z = [u, T , v, T ] T ∈ D(A). Thus, A is dissipative. Let the norm induced by (43) square be the Lyapunov function W , that is, W = ⟨z, z⟩ E = ‖z‖ 2 E . Since the time derivative of W along the solution of the system is given by dW ∕dt = 2⟨Az, z⟩ H , the operator is dissipative means dW ∕dt ≤ 0.
Next, we show 0 ∈ (A), where (A) is the resolvent set of the operator A. For this proof, we have to show that −A −1 exists and is bounded. For any given = [ 1 , 2 , … , 14 ] T ∈ H , we find = [ 1 , 2 , … , 14 ] T ∈ D(A) so that A = − . Directly solving A = − , we obtain The coefficient c i , i = 0, … , 3, can be determined by boundary conditions, and we obtain the matrix form relation to (49), we see that there exists K > 0 such that Thus, we obtain 0 ∈ (A). As the operator A is dissipative and 0 ∈ (A), the operator A generates a C 0 -semi-group of contractions from the Lumer-Phillips theorem [35]. Further, we can show the compactness of A −1 . In fact, we obtain the following inequality by evaluating the upper bound of in Equation (49): where K ′ is a positive constant. By applying Sobolev embedding theorem [36], the inverse operator of A is compact on H . Therefore LaSalle's invariance principle is applicable to the system (46) [37]. □

Asymptotic stability
From Lemma 1, we found that all solutions of the system (46) asymptotically converge to the maximal invariant subset of the following set  : Now, we show that the set  contains only zero solution, that is the desired states.

Lemma 2. The set  contains only desired states.
Proof. Under the condition dW ∕dt = 0, the closed-loop system (38)-(40) becomes following: In the following, a solution is obtained by using the separation of variables. We assume solution can be expressed as w(x, t ) = X (x)T (t ), where X (x) and T (t ) denote functions of space and time, respectively. Now, we consider the following two cases; T (t ) ≠ 0 and T (t ) = 0.

T (t ) ≠ 0
Equation (54) can be written as follows: where = const., and we obtain X ′′′′ (x) + X (x) = 0 for X (x) andT (t ) = 2 T (t ), 2 = EI for T (t ). Now, we consider three cases depends on values of . (a) = − 4 < 0; general solution of T (t ) is and this contradicts that the energy is finite.

Theorem 1. The system (46) is asymptotically stable.
Proof. Let (t ) be a C 0 -semi-group of contractions generated by operator A (44). Lemma 1 means that the system (46) has a unique solution z(t ) = (t )z 0 , for the initial value z 0 ∈ D(A). For such a system, we can show the asymptotic stability using the extended LaSalle's invariance principle [37,Th.4.3.4]. According to this principle, the solution of (46) converges to the maximal invariant set of  , provided that the solution trajectories are pre-compact. Here, we can see that the operator ( I − A) −1 compact for ∈ (A) from the fact that A −1 is compact and [38,Th.6.29]. Thus, the solution trajectories of (46) are pre-compact [39,Th.5.2]. Further, the maximal invariant set of  contains only desired states from Lemma 2. Therefore, the system at the neighbourhood of the desired configuration is asymptotically stable. □

Simulation setup
The effectiveness of the proposed controller is investigated by numerical simulations. To obtain the finite-dimensional model of the flexible structure, we used a constrained mode method [40] and considered the first seven vibration modes. For any point on the beam, the elastic deformation w(x, t ) is approximated as where i (x) is the ith mode function of the beam, i is the eigenvalue correspond to the "clamped-free" boundary conditions, p i (t ) is ith generalised coordinate representing the contribution of the ith mode shape. Substituting Equation (58) into Equation (16), we obtain the dynamic equation of p i (t ) by using the orthogonal characteristics of i (x). For more details of the finitedimensional model of a flexible beam, see [41,42]. The nonlinear equations of motion have been integrated via a fourth-order Runge-Kutta method. In addition, parameters of the flexible beam and robots are shown in Table 2. We set sampling time In the simulation, we carried out three simulations. The first (Simulation 1) is the step response, the second (Simulation 2) is the disturbance response, and the third (Simulation 3) is the step response when the thrusters can generate on/off thrust only. For the step responses (Simulation 1 and 3), we set the desired configuration as follows: r Xd = r Yd = 10, (25). On the other hand, for the disturbance response (Simulation 2), we applied the pulse disturbances −10 (N) to F ix (t ), F iy (t ), and i (t ), i = 1, 2 from 50 (s) to 51 (s) as a disturbance caused by the thruster. In the disturbance response, the desired configuration was set to be the same as the initial conditions. In each case, the initial condition is [r X (0), r Y (0), (0), q 11 (0), q 21 (0), w(x, 0)]= 0 and [̇r X (0),̇r Y (0),̇(0),̇q 11 (0),̇q 21 (0),ẇ(x, 0)]= 0. The control parameter is k X = k Y = 3∕2, k = 30, k 11 = k 21 = 2∕5, k 1 = k 2 = 20, k 3 = 2 × 10 3 , k 4 = k 5 = 10, k 6 = 2, k 7 = 5. Grasping angle q ib , i = 1, 2, is ∕2 (rad). In the step response, we also compare the input values with and without using the null space. Weighting matrix is set as W 1 = diag[1, 1, 1∕2, 1∕2, 1, 1, 1∕2, 1∕2] which intended to reduce the input of the robot1 to the robot2.

Simulation 1
First, we investigated the step response. Figure 2 shows  Figure 3 shows the step response of w(x, t ) by the proposed controller, and the dotted line denotes the initial value of the transverse displacement. From these figures, state variables were converged to desired values, and the transverse displacement of the beam was effectively suppressed in the proposed method. Thus, the proposed controller could realise the control objective. Next, we compared the input values with and without using actuator redundancy discussed in sub-section 2.2. The input with using actuator redundancy is Equation (24), and without using actuator redundancy is Equation (21) with = 0.  show the responses of the proposed controller using actuator redundancy and the proposed controller without using actuator redundancy, respectively. From the results, the norm of robot 1's inputs is effectively reduced by using actuator redundancy. However, as the trade-off, the input norm of robot 2 increased. In the simulation, the weighting matrix was set as W 1 = diag[1, 1, 1∕2, 1∕2, 1, 1, 1∕2, 1∕2], and thus the input of the robot 1 was reduced compared to the input of the robot 2. This is useful when we need to save the robot 1's thruster fuel and battery.

Simulation 2
Next, we investigated the disturbance response. Here, we compared the response of the system with the proposed controller and without the controller. For the system without controller, we set u = 0. Figure 7 shows the responses of state variables of the robots' state, and Figure 8 shows the time response of w(L, t ). In both figures, a solid line and a dashed line show the responses of the system with the proposed controller and without the controller, respectively. In Figure 8, the system's response with the proposed controller was quickly converged to the desired value, but the system's response without controller left some vibrations. In addition, from Figure 7, the system's responses without controller did not converge to the desired values, and their values largely deviated from the desired values. On the other hand, the system's responses with the proposed controller were converged to the desired values. From these points of view, we found that

Simulation 3
To investigate the more practical applicability of the proposed controller, we examined the step response when the thrusters can generate on/off thrust only. Here, we assumed the following for the thrusters' output: the thrusters generate on/off thrust; and the output has the maximum value of T out . Thruster on/off switching time is infinitesimally small. In addition, due to the combination of multiple thrusters, the translation forces by thrusters on Robot i' base, F ix (t ) and F iy (t ), take a value of 0 or +T out , or −T out . In particular, if the current state is far from the desired value, a large output thruster (maximum output T out = 20 N) is used, and if the current state is close to the desired value, a small output thruster (maximum output T out = 1 N) is used. The change of the thruster used is as follows. If |r X (t ) − r Xd | ≤ 2 and |r Y (t ) − r Yd | ≤ 2 are satisfied, a small output thruster is used. If this condition is not satisfied, a large output thruster is used. A feedback modulator [43,44], was used to convert the control input Equations (24), (36), which is a continuous signal, into an on/off signal. The configuration of the feedback modulator is shown in Figure 9. In this figure, S/H denotes the zeroorder holder with the sample period h = 0.5 (s), and Quantiser denotes the quantiser. In this paper, to generate an on/off output, this quantiser generates a value of 0 or +T out , or −T out . In addition, u means the continuous control input, u c is the input Transverse displacement w(x, t ) by proposed controller with on/off thrusters of S/H, and u Q is the on/off output. The system's response becomes unstable just by converting the continuous signal into the on/off signal using S/H and the quantiser. The feedback modulator feeds back the error due to S/H and the quantiser to the input signal u. As a result, we can convert the continuous control input to an on/off input signal while maintaining the closed-loop system's stability. Here, Q(s) = 1 − ( s∕( s + 1)) 2 is a design parameter, and the stability can be maintained by choosing > h [43]. Here, we set = 0.501. Figures 10-12 show the step response of state variables related to the robots' state, w(x, t ), and time response of input  Figure 12, we can see that translation forces by thrusters on Robot i's base, F ix (t ) and F iy (t ), take a value of 0 or +T out , or −T out . Particularly, the large thruster was used for up to 20 seconds, then the usage was once changed to the smaller thruster, and the larger thruster was used again. After 50 seconds, the small thruster was used. On the other hand, by comparing Figures 2 and 3 with Figures 10 and 11, we found that the system response slows down to the desired value due to on/off thrusters. However, all responses converge to the desired value without divergence, and we could confirm the practical applicability of the proposed controller.
The impulse modulators commonly used have Schmitt trigger, pulse-width modulator (PWM), pulse-width pulsefrequency (PWPF) modulator [45,46]. However, these PWM will not achieve the desired performance unless its pulse period can be made sufficiently shorter than the control period. On the other hand, the feedback modulator feeds back the error due to S/H and the quantiser, and the modulator suppress their influences. Thus, the feedback modulator works well with slow pulse periods, unlike the PWM case. Because of this convenience, we used a feedback modulator in this paper. A conventional modulator, PWM, and a feedback modulator are compared in [43], and see it for details.

CONCLUSIONS
In this paper, we considered a cooperative transportation control problem of flexible Euler-Bernoulli beam by two one-link rigid planar space robots. To solve the problem, we proposed a boundary controller based on the hybrid PDE-ODE model. The proposed controller has a simple structure and lets the closed-loop system become Lyapunov stable. Further, because of the original system is too complicated to prove global asymptotic stability, we proved that the proposed controller could ensure asymptotic stability at the neighbourhood of the desired configuration. In addition, we designed the input by considering each actuator's performance because of the actuator redundancy of the system. Finally, in numerical simulations, we confirmed that the proposed controller works well for the step responses and the disturbance responses.
In this paper, we consider the point to point transportation. To make the controller more practical, trajectory tracking control will be desirable as future research. In addition, the theoretical analysis of the disturbance response is essential, and we plan to address it in the future.