A new computational approach for optimal control of switched systems

The combination of the time-scaling transformation and control parameterization has proven to be an eﬀective approach in addressing optimal control problems involving switching systems with predeﬁned subsystem sequences. However, this approach has certain limitations. First, the number of control switchings is required to be no less than the number of subsystem switchings. Second, the switching of the subsystem must be accompanied by the switching of the control. Third, this scheme introduces many hyperparameters, leading to combinatorial explosion. To address these drawbacks, we introduce a novel computational approach such that the control switching can be independent of subsystem switching. The superiority of this novel approach can be clearly observed from the solutions obtained using the proposed method for solving two illustrative examples


Introduction
A switching system is a dynamic system composed of multiple modes (often called subsystems).The system transitions between these modes based on defined rules or conditions.More specifically, a transition between subsystems can be triggered by predefined events, specific time instances, or conditions tied to the system state.Due to their inherent adaptability, such systems can dynamically shape their behavior and control strategies to adeptly respond to changes in environmental conditions or system requirements.Switching systems are encountered in a variety of areas such as robots [1,2], power systems [3], and communication networks [4,5].The success of these applications demonstrates the versatility and importance of switching systems.The challenge in solving the optimal control problem of switching systems (SOCP) lies in intelligently determining the control inputs and switching laws.The switching law serves as a blueprint outlining the detailed sequence of active subsystems and their corresponding switching times.This sequence encapsulates the precise moments when the system transitions smoothly from one subsystem to another, coordinating the interaction of modes to optimize overall performance.
Early theoretical advances to address this problem can be found in [6] and [7], where dynamic programming and the maximum principle were adopted, respectively.However, they do not provide effective numerical methods to obtain optimal solutions.In [8], Xu and Antsaklis introduced a bilevel optimization approach to address SOCPs.This approach breaks down the problem into two levels and solves them iteratively.The lower level focuses on fixing the switching sequence and optimizing the control inputs and switching times to minimize the cost function.On the other hand, the upper layer optimizes the switching sequence and performs the lower layer optimization for each switching sequence.This is much more complex.In [9][10][11][12][13], it was generally assumed that the switching sequence is predetermined, and the focus is on optimizing the lower level.This consideration applies to scenarios where the order of switching sequences is predetermined and only lower-level issues are performed.See, for example, [8].In addition, the research results of the lower level will lay the foundation for the research of the upper level.Therefore, this paper also studies the switching optimal control problem of fixed switching sequence.
Most existing strategies for addressing the lower-level problem resort to discretizing the control function and/or state variable.Then, a search algorithm is applied to the resulting discrete model to find the optimal solution [14].However, these discretization methods can lead to combinatorial explosion.To address this issue, a novel method is proposed in [8] to decompose the lower-level problem into two subproblems.The first subproblem focuses on minimizing the cost function with fixed switching instants, where decision variables are limited to control inputs.This process converts the problem into a nonlinear optimization problem.For the second subproblem, it aims to minimize the cost function with respect to switching instants.For this, the gradient of the cost function with respect to the switching times is derived.This approach has been improved in [15][16][17].
In [18], a novel technique is introduced to address the lower-level problem without the need for decomposition.This method works by first approximating the control function using a piecewise constant function where the control heights and the control switching times become the new optimization variables.Then, with the aid of the time-scaling transformation [19,20], the switching times of the control and subsystems are mapped into fixed knots within a newly defined time horizon.This process leads to the approximation of the original optimal control problem as a nonlinear optimization problem.Since then, this technique has been widely adopted for solving optimal control problems involving switching systems over the past decade [21][22][23].However, it is worth noting that this scheme requires that (i) the number of control switchings must be no less than the number of switchings between subsystems, and (ii) the switching of the subsystem must be accompanied by the switching of the control.These predetermined requirements may not be strictly adhered to for many practical problems.This is because the control switching times and the subsystem switching times are usually independent of each other.
In [24,25], within the computational paradigm of control parameterization, an effective strategy called Sequential Adaptive Switching Time Optimization (SASTO) is proposed to solve standard optimal control problems.The SASTO approach is applicable to scenarios where the optimization is to be carried out on multiple sets of mutually independent switching times.It is achieved through the systematic introduction of new time domains.Motivated by the success of this approach, our aim is to devise a new method that efficiently addresses the limitations faced when solving switching optimal control problems.
The remainder of the paper is structured as follows.The problem formulation is presented in Sect. 2. Moving on to Sect. 3, we delve into the shortcomings of employing traditional time scale transformations and devise an efficient approach to address this challenge.A new computational method is constructed to solve the resulting transformed problem.In Sect.4, the gradient formulas of the cost/constraint functions are given.In Sect.5, two numerical examples are solved.Finally, concluding remarks are provided in Sect.6.

Problem formulation
Consider the following nonlinear switched system with N subsystems: where x(t) ∈ R n and u(t) ∈ R m are the state and control vectors of the system at time t, respectively; x 0 ∈ R n is the given initial state; f j : R n × R m → R n , j = 1, . . ., N , is the given continuously differentiable function representing the jth subsystem; t 0 , t 1 , . . ., t N are the switching instants of subsystems where t 0 = 0, t N = T > 0 are the initial and terminal times, respectively.Let τ = [t 0 , t 1 , . . . ,t N ] be such that t j-1 < t j , j = 1, . . . ,N .Let T be the set consisting of all such τ .Let u : [0, T] → R m be a Borel measurable function satisfying u(t) ∈ U for almost all t ∈ [0, T], where U is a compact and convex set in R m .Such a u is termed an admissible control.Denote the set of all such admissible controls as U .The formal definition of the optimal control problem for the switched system is as follows: Problem (P1) Given the switched system (1)-( 2), find u ∈ U and τ ∈ T such that the cost function Here, N b and N d represent the counts of inequality and equality constraints, respectively; We presume that the conditions outlined in [20] are satisfied throughout this paper: A1.Function f j , j = 1, . . ., N , is twice continuously differentiable.A2.For j = 1, . . ., N , there is a constant L (L > 0) such that A3. Functions l : R n → R, l = 0, 1, . . . ,N b + N d , and L l : R n × R m → R, l = 0, 1, . . ., N b + N d , are continuously differentiable with respect to their arguments.

Method analysis
Over the past few decades, the utilization of the time-scaling transformation and control parameterization method has gained popularity as a numerical method for solving SOCPs [18,21,26].This method subdivides the time interval [t j-1 , t j ] into p j subintervals and denotes the p j + 1 nodes as τ j 0 , τ j 1 , . . ., τ j p j , j = 1, . . ., N , where t 0 = 0, t N = T and as shown in Fig. 1.By applying the control parameterization, u(t) is approximated by a piecewise constant function on the time domain [0, T] as follows: where ρ j,k is the approximate value of the control u(t) on the kth interval [t j-1 , t j ].Then, the time-scaling transformation operates by mapping these variable knots {t 0 , τ 1 1 , . . ., Clearly, there are many hyperparameters p 1 , . . ., p N , which can easily lead to combinatoric explosion.For example, when p = 5 and p N j=1 p j , there are 21 possible combinations of partitions (p 1 , p 2 , p 3 ).To avoid this situation, most existing literature such as [18,21,22,27,28] only consider the case of p 1 = p 2 = p 3 in their numerical experiments so that p must be a multiple of N .
Additionally, from Fig. 1 it can be seen that the switching instants t 0 , t 1 , . . ., t N coincide with the partition knots τ 1 0 , τ 2 0 , . . ., τ N N , so p N j=1 p j ≥ N , which indicates (i) the number of control switchings must be no less than the number of switchings between subsystems; and (ii) the switching of the subsystem must be accompanied by the switching of the control.This is clearly unreasonable because the switching times t 0 , t 1 , . . ., t N and the control switching points are usually mutually independent.
In [24,25], it is pointed out that the time-scaling transformation is an order-preserving mapping.Therefore, the application of this method must meet the essential requirements of a priori knowledge of the chronological order of these time instants to be optimized.
Nevertheless, due to the independence of these two types of switching instants, the temporal sequence between the control switching points and the switching times of the subsystems cannot be known in advance.Thus, using this method to solve Problem (P1) must require the switching of the subsystem be accompanied by the switching of control, i.e., t j = τ j+1 0 , j = 1, . . ., N .In other words, the prior knowledge of the chronological order of these time instants {t 0 , τ 1  1 , . . ., This clearly violates the independent nature between these two types of switching instants.
To overcome these shortcomings, we will design a new method to solve Problem (P1), which enables adaptive optimization of the selection of unknown switching time sequences as in [24].

A new numerical solution procedure
For Problem (P1), the switching instant vector τ and the control function u(t) need to be optimized jointly.Initially, we employ the time-scaling transformation to the subsystem which involves τ , while keeping u(t) temporarily unchanged.This process leads to the formulation of a new optimal control problem.Following this, we implement the control parameterization for u(t) within the newly defined time domain and subsequently utilize the time-scaling transformation to handle the control switching times, thereby introducing an additional new time domain.The detailed procedure is given below.
Step 1. Time-scaling transformation for switching time vector τ For each τ = [t 1 , . . . ,t N ] ∈ T , define θ = [θ 1 , . . ., θ N ] ∈ R N , where θ j = t jt j-1 ≥ 0, j = 1, . . ., N , and θ 1 + • • •+ θ N = T. Define as the set comprising all such vectors θ ∈ R N .By introducing a new time variable s, we establish the definition for the time-scaling function Here, • represents the floor function.The choice of θ N+1 is arbitrary since its value does not impact the result when s = N .Under the effect of time-scaling transformation, t = t j are transformed into fixed points s = j, j = 1, . . ., N .Figure 2 shows the transformation of the variable switching times to the fixed points on the time domain s through the time-scaling function μ 1 (s | θ ) when N = 4. Accordingly, u(t) is transformed into u(μ 1 (s)) w(s), which is called the transitional control function of u(t).The switching system defined on the new time domain [0, N] is ẏ(s) = θ j f j y(s), w(s) , s ∈ [j -1, j), j = 1, . . ., N, ( 4 ) where y(s) x(μ 1 (s)) = x(t) is the new state vector of the time variable s.Step 2. Control parameterization applied to w(s) We divide the time interval [0, N] into p subintervals with p + 1 partition points denoted by s 0 , s 1 , . . ., s p , satisfying 0 = s 0 < s 1 < • • • < s p = N as shown in Fig. 3.Note that we introduce only one hyperparameter p compared with the method in [18] (see Fig. 1).Define as the set comprising all such vectors ξ = [s 0 , s 1 , . . ., s p ] ∈ R p , where s k , k = 1, . . ., p, are the switching times of the transitional control function w(s).Accordingly, w(s) is approximated by where is the approximate value of the transitional control w(s) = [w 1 (s), . . ., w m (s)] on the kth subinterval.Define as the set comprising all such vectors ρ = [ρ 1 , . . ., ρ p ] and let χ [s k-1 ,s k ) (s) be defined by The switching system (4)-( 5) defined on the interval [0, N] becomes ẏ(s) = θ j f j y(s), Step 3. Time-scaling transformation applied to the control switching times in ξ For each ξ = [s 0 , s 1 , . . ., The above transformation maps variable times s = s k to fixed times v = k, k = 1, . . ., p. Also y(s) = y(μ 2 (v)) z(v), which is the new state vector of the time variable v. Accordingly, the switching system (7)-(8) becomes Denote z(•|θ , σ , ρ) as the solution of ( 9)- (10).Based on Assumptions A1-A2, we can show that the system possesses a unique solution for every permissible (θ , σ , ρ) within the set ( , , ) [20].Problem (P1) is approximated as the following nonlinear programming problem.

The switching times analysis
After making the above approximation to Problem (P1), the original switching time vector τ is transformed into the duration θ between two successive subsystem switching points, while the switching times for the control function u(t) are transformed into the control height ρ and the duration σ between two successive control switching knots.Here, ρ, σ , and θ are new decision vectors to be optimized.Once the optimal solution (ρ * , σ * , θ * ) is obtained by solving Problem (P2), we can derive the switching times of the subsystem and the control on the original time interval [0, T].
To be specific, the optimal switching times, t * j , j = 1, . . ., N , of the subsystem can be readily calculated by using the time-scaling function μ 1 (s | θ * ) and setting s = 1, . . ., N .However, determining the switching times of u(t) necessitates the following procedure: . ., τ * p ] are selected adaptively, being carried out under the condition that the subsystem switching times and the control switching times are mutually independent, as illustrated in Fig. 4.Moreover, it is worth noting that compared to Fig. 1, there is only one hyperparameter p.

Figure 4
The optimal switching times of the subsystems and the control

Gradient computation
Problem (P2) can be addressed through the application of gradient-based nonlinear optimization algorithms.To accomplish this, it is necessary to calculate the gradients of both the cost and constraint functions, where the used gradient computational method is the variational method.The theorem below provides the partial derivatives of the state at θ .
Proof The proof closely resembles that presented in Theorem 1 of [29] and is consequently omitted.

Theorem 2
For each l = 0, 1, . . ., N b + N d , the gradient of g l in terms of θ , are given by Proof The proof is derived by using the chain rule to the cost/constraint functionals in conjunction with Theorem 1.
Likewise, the derivatives of the cost/constraint functions in terms of ρ can be computed. where Utilizing the chain rule, the theorem below provides the gradients of g l (θ , σ , ρ) at ρ.

Theorem 4
For each l = 0, 1, . . ., N b + N d , the gradient of g l in terms of ρ, are given by Regarding the gradients of the cost/constraint functions in terms of σ , the derivation is more intricate, as discussed in [24].To begin, we provide the inverse function of where k = 2, . . ., p.The switching times of subsystems in the time domain v, denoted as μ -1 2 (s | σ ), are dependent on σ .At each junction point μ -1 2 (h) where h = 1, . . ., N -1, coinciding with any v = 1, . . ., p -1, the inverse function is nondifferentiable.

Case 1. When the point v h aligns with any
Case 2. We have ∂σ when v h does not align with any v = 1, . . ., p -1.

Theorem 5
For each (θ , σ , ρ), where (•|θ , σ , ρ) is the solution of the following variational system on v ∈ [k -1, k), k = 1, . . ., p: Furthermore, for each l = 0, 1, . . ., N b + N d , the partial derivative of ḡl in terms of σ is given by The definitions of H - h,L l and H + h,L l are analogous to those of H - h,f and H + h,f .

Illustrative example
To showcase the efficacy of the proposed approach, we examine two representative examples.We employ the Fmincon-SQP optimization solver in MATLAB R2022a to address our problem.The ordinary differential equations (ODEs) are solved using the fourth-order Runge-Kutta method with a step size of 10 -2 .The calculation of cost and constraint integrals is performed using Simpson's rule, with a step size of 10 -2 as well.All numerical experiments are conducted on a computer with a CPU and RAM of 3.60 GHz and 16 GB, respectively.
Example 5.1 Consider a nonlinear switching system, as detailed in [8], composed of three subsystems described as follows: subsystem 2: subsystem 3: Let t 0 = 0 and t f = 3, and assume that the system switches at t = t 1 from subsystem 1 to subsystem 2 and at t = t 2 from subsystem 2 to subsystem 3 (0 < t 1 < t 2 < 3).Our objective is to determine the optimal switching instants t 1 , t 2 and the optimal input u to minimize the cost functional Here, x(0) = [1,1] .First, we apply the new approach to solve this problem using 3, 4, and 5 partitions for the control function within the time interval of [0,3].For comparison, we also consider using the conventional time-scaling transformation to solve this problem.As stated in Sect.3, this method requires that p 1 = p 2 = p 3 (p 3 j=1 p j must be a multiple of N ).However, a crucial distinction between them is that the time-scaling method requires that the system switching times and discrete time knots of the control space coincide, resulting in the computational combinatoric explosion.On the other hand, the proposed method can effectively avoid the above deficiency by allowing the system switching to be independent of the control switching.To continue, when applying the time-scaling method, we select p = 3, 6, 9 with equal partition numbers (p 1 = p 2 = p 3 ).The corresponding optimal results obtained using these two methods are listed in Table 1.The optimal control are shown in Figs.5-7, where the two red vertical lines represent t = t 1 and t = t 2 , respectively; t 1 and t 2 are the optimal switching times of the subsystem.
Table 1 reveals that the optimal cost, 3.6809, obtained by doubling p from 3 to 9 using the time-scaling method is inferior to the result obtained by using the proposed method  for p = 5, which indicates the effectiveness of our method in eliminating unnecessary partitions.Figure 5 shows that the time-scaling method necessitates the optimal control switches at t 1 and t 2 , while our method has no such restriction.Furthermore, in Fig. 7, the optimal number of switches of the control u on each subinterval [t i-1 , t i ], i = 1, 2, 3, is given in advance and is equal.However, in Fig. 6, the optimal number of switches of the control u on each subinterval is determined by the adaptive optimization.This demonstrates our method can achieve the adaptive selection of the number of control segments within each subinterval [t i-1 , t i ], i = 1, 2, 3, whereas the time-scaling method lacks this flexibility.
Example 5.2 In this problem [18], a nonlinear switching system with three subsystems is defined over the time domain [0, 3]: subsystem 1: subsystem 2: The initial conditions are The system undergoes a switch at t = t 1 from subsystem 1 to subsystem 2 and at t = t 2 from subsystem 2 to subsystem 3 (0 < t 1 < t 2 < 3).We aim to determine the optimal switching instants t 1 , t 2 , and the optimal input u that minimize the cost function Note that the cost function contains the state x 1 (t 1 ) at the switching time t = t 1 .
We use the proposed method and the time-scaling technique to solve this problem.The optimal results are listed in Table 2.The optimal controls are depicted in Figs.8-9, with the two red vertical lines denoting t = t 1 and t = t 2 correspondingly.Here, t 1 and t 2 are the optimal switching times of the subsystem.
It can be seen from Table 2 that when p = 6 and p = 9, the optimal costs obtained by the proposed method both are obviously better than the optimal costs obtained by the time-scaling method.Figure 8 when compared with Fig. 9 also shows that our method can avoid the limitation that the switching of the subsystem must be accompanied by the switching of the control, and allows for the adaptive determination of the number of control segments within each subinterval.

Conclusion
We considered a type of optimal control problems involving the switched system with a prespecified sequence of subsystems.We then highlighted the limitations of the conventional time-scaling transformation used in solving this type of optimal control problems.Subsequently, we propose a novel method to overcome the existing shortcomings of the conventional time scaling transformation, where the control switching times and system switching moments are sequentially selected independently.Notably, in the proposed method, we introduced only one hyperparameter, reducing the need for extensive parameter tuning, thus improving overall efficiency.Two numerical examples are solved and the results obtained clearly illustrate the effectiveness of the new approach.In future studies, the proposed method can also be considered to solve problems involving the switching time-delay system and the switching time-dependent state system.

Figure 1
Figure 1 Illustration of switching times for u(t)

Figure 5 1 Figure 6 1 Figure 7
Figure 5 Optimal control using these two method when p = 3 for Example 5.1

Figure 8 2 Figure 9
Figure 8Optimal control using the proposed method when p = 6 for Example 5.2

Table 1
Optimal results and computational times using the two techniques for Example 5.1

Table 2
Optimal results and computational times using the two techniques for Example 5.2