A global optimization approach for sliding mode tuning and existence maps generation

In this work, global optimization techniques based on interval arithmetic are proposed to analyze and synthesize sliding mode (SM) controllers. The proposed methodology allows generating a series of maps, called subpavings, which put in evidence the required relationships among tuning parameters, disturbances and control amplitude to fulfill the sufficient condition of SM in a guaranteed way. The a priori knowledge of the control power necessary to guarantee SM behavior of nonlinear systems despite parameter uncertainties and external disturbances is an advantage of this proposal compared to traditional tuning methods, which usually fall into over-sizing solutions. Although the methodology is developed in the context of conventional first-order SM controllers, it could be extended to any other SM design approach.


Introduction
The design of controllers based on sliding modes (SM) is characterized by its applicability to nonlinear systems and by its robustness. Extensive studies have been carried out in this area based on classical methods like state-feedback sliding surface design, optimization-based designs or even SM methods that combine adaptive control features [1][2][3].
Conventional designs are performed in two stages. The first stage is devoted to find a sliding manifold on which the desired closed-loop dynamics is achieved. The second stage consists in designing a discontinuous action and a switching logic that enforce the system to reach the prescribed manifold and to slide on it from then on [4][5][6]. The challenge in this type of design is to find a manifold and a control action that satisfy the necessary and sufficient condition for SM existence, at least in the desired operating range of the system. The traditional design consists in delimiting the parameters of the system by its extreme values, and then to obtain the constant control action necessary for SM operation which is finally tested through simulation [1,7,8]. However, this non-optimal design may not guarantee SM for some possible states of the system, or even worse, it could lead to extremely conservative solutions thus degrading the controller performance (according to, for instance, power consumption or chattering).
To overcome these problems, adaptive and optimization techniques have been exploited. Adaptive methods adjust online the coefficients of the chosen surface and the control action with the aim of reducing chattering and unmatched disturbances effects. Azimi and Shahravi [9] and Wang et al. [10] use a time-varying commutation law, where its coefficients are adapted online and the exact knowledge of the system is not a priori necessary. Utkin et al. [11] proposes a similar method seeking to reduce the value of discontinuous action in order to mitigate chattering effects. Wan et al. [12] combine SM design with a neural network that adapts the control signal according to its distance to the desired surface. In all of these works, although stability and robustness criteria are given, the robustness margin is not generally analysed, i.e. up to which excursion in system parameters the expected behavior will remain. Optimization methods, on the other hand, seek to optimize a functional via quadratic optimization or linear matrix inequalities, looking for the optimal sliding manifold. In Shtessel et al. [8] the optimal tuning of the SM manifold is transformed into a linear quadratic regulation (LQR) problem, while the design of the discontinuous control action is addressed separately. In contrast, Edwards [13] provides an optimal design, via LMIs, that considers both the SM dynamics and the control costs required to maintain sliding. Also, the work of Salamci et al. [14] seeks to solve the problem of SM optimal controllers in nonlinear systems using a recursive approximation with linear models together with an optimization criterion for the selection of the control to apply. Despite their attractive features, these approaches apply to linear systems or require linearization in order to perform the optimization.
In this paper, an alternative design methodology is proposed, which essentially inverts the order of the design steps. That is, starting from the system dynamics, bounded uncertainties and disturbances, and control power, an optimal SM dynamics, i.e. an optimal sliding manifold, is derived. To this end, the proposal integrates optimization techniques based on interval arithmetic with classical SM tools. Differing from other proposals, global optimization techniques are used here so that the approach can be applied with guarantees to nonlinear systems. In particular, Interval Branch and Bound (IBBA) and Set Inversion Via Interval Analysis (SIVIA) algorithms are combined to find an optimal controller for a particular criterion, considering the possible uncertainties and disturbances of the nonlinear system. Also, a series of maps called subpavings are obtained that provide the state-space subsets where the designed control will have a guaranteed behavior. These maps serve as SM design tool even for adaptive approaches. Summarizing, the main contributions of the paper are: -An SM optimal tuning controller methodology that uses a customizable optimization function and can deal with nonlinear dynamics subjected to parameter uncertainties and external disturbances. -A series of maps that determine areas of guaranteed SM operation and show sensitivity to parameter variations, thus serving as analysis and design tools for SM controlled systems.
As study case, the proposed tuning method is applied to the heave direction of the autonomous underwater vehicle (AUV) Ciscrea (see Fig. 1). The environment where AUVs work are unstructured, partially unknown and highly perturbed. Furthermore, AUVs are very sensitive to model-based controller designs because their models include uncertain hydrodynamic parameters [15][16][17]. All these challenges justify robust control approaches like SM techniques with particular care about system's uncertainties and disturbances. Several papers in the bibliography have applied SM control to deal with the challenges of the application, particularly using adaptive algorithms [3,18,19]. The use of interval arithmetic in the context of AUVs has been more relegated to state estimation [20,21]. In the context of SM, its use has been as a complement to get on-line controllers with continuous adaptation of the command signal [22,23], or to add robustness to the controller design [24]. Nevertheless, global optimization based on interval arithmetic can also be used to synthesize controllers from different criteria or control methods, such as H ∞ criteria [25] or Quantitative Feedback Theory [26]. One of the objectives of this paper is to extend the global optimization approach to sliding mode control.
The article is structured as follows: Sect. 2 establishes the principles of SM and the arithmetic of intervals, as well as the statement of the optimization problem. Section 3 formulates SM analysis and synthesis problems as optimization ones. Then, Sect. 4 applies the optimization approach to the SM control of the AUV Ciscrea (see Fig. 1). Finally, some concluding remarks are given.

Theoretical framework
In this section an introduction to SM, global optimization and interval analysis is performed. The concepts given here will provide the basis for understanding the Sect. 3 "Proposal: GO as a SM tuning tool".

SM control essentials
In order to design an SM controller, a manifold in the state space must be selected according to some performance criterion. This manifold, which must be in accordance with the physical constraints of the system, governs the systems dynamics during the sliding motion. In addition to the manifold selection, a discontinuous law must be designed to enforce the state trajectories towards the manifold from both sides.
Consider a nonlinear perturbed system of the form: where x ∈ R n is the system state, u ∈ R is the control signal, y ∈ R is the system output, f (x) and g(x) are the drift and control vector fields in R n , d(x, t) ∈ R n × R + is the disturbance vector encompassing uncertain parameters and external perturbations, and finally h(x) is a scalar function in R n . The switching policy can be defined as: in accordance with the sign of the auxiliary output σ (x), which can be interpreted as the system output that vanishes on the desired manifold The sliding surface and switching law must be designed to fulfill the so-called reaching condition that is sufficient and necessary for SM existence [1]: For this condition to locally hold on both sides of S,σ (x) must depend on the discontinuous signal u. This necessary condition constrains the family of sliding surface candidates. Additionally, u − and u + must be properly designed. After reaching the sliding surface S, provided (4) holds, a very high frequency (ideally infinite) commutation is established, forcing the system to slide on it.
There exists a smooth control equivalent to the discontinuous law in the sense that, once on the surface S, the system state moves on S in the same way as it does with the discontinuous control [1]. This so-called equivalent control u eq (x, t) can be obtained from: where the operator L f h(x) : R n → R (directional or Lie derivative) denotes the derivative of a scalar field h(x) : (5), the smooth control law that makes S an invariant subset is given by: Taking this into account, the necessary and sufficient condition for SM existence (4) can be rewritten in terms of u eq . From (6), L g σ = 0 establishes the necessary condition for the existence of u eq and SM. Moreover, a sufficient condition for the local existence of the SM over S can be derived from (4) and (5), considering (without loss of generality) u + > u − : Obviously, this equivalent control is fictitious and in general cannot be constructed because of its dependence on the uncertain d(x, t).
With regards to the disturbances, d(x, t) can always be decomposed into two terms, one belonging to the tangent space of S and the other aligned with g(x). Perturbations that are tangent to S alter the SM dynamics but not its existence condition. On the other hand, the SM dynamics is completely invariant to perturbations aligned with g(x), which affect the equivalent control and therefore the existence condition. In fact, for a matched perturbation d( is the equivalent control of the unperturbed system. Furthermore, if uncertainty in the control vector field is also considered, i.e. g( is the control vector field of the unperturbed system, then u eq is related with the unperturbed u * eq as follows (see [1] and references therein): Whenever the necessary and sufficient condition for SM existence (7) holds in a given region of the state space for any bounded disturbance, the controller will be said robust. If, in addition, there is no disturbance on the tangent subspace of the sliding surface, i.e. if the disturbance is aligned with g(x), the SM will be said strongly invariant.

Global optimization (GO)
The optimization goal is to achieve an optimal design of the sliding surface for a given criterion. This requires checking the necessary and sufficient condition of SM existence for all admissible excursions of state variables, system parameters and external disturbances. With this aim, consider a continuous constrained optimization problem, where is the function that defines the subset where the solution is searched.
The solution to this problem (the minimizer), is expressed as k * . It is the point where m is minimum over the set When m and c are not convex functions, local optimization techniques do not guarantee convergence to the global solution k * . In contrast, global optimization methods converge to the global minimum and supply an enclosure [m * , m * ] of m * .
In the following subsection, the global optimization algorithm of Branch and Bound based on interval arithmetic is revisited [27].

Interval arithmetic
The proposal of this work is based on the Interval Branch and Bound Algorithm (IBBA) and the Set Inversion Via Interval Analysis (SIVIA) algorithms [28]. The concepts of real intervals and inclusion functions are introduced here as they will be used to represent the bounded uncertainties and sliding functions in the next section. This way, the IBBA algorithm will be exploited as a robust SM tuning tool.
is an n-dimensional interval vector belonging to the space IR n [29].
Considering the problem (9), if inclusion functions of m and c exist, k * can be searched in K ⊂ IR n . Furthermore, a guaranteed lower bound m and an upper bound m of m * can be achieved using the IBBA. To compute this, IBBA repeatedly bisects K in smaller boxes [k i ] and discards them when it is proven that or if a feasible pointk has been found (through the testing of random points in each box) such that any points in [k i ] can provide a better feasible solution, The IBBA stops when the distance between m and m reaches the desired precision , with Figure 2 illustrates the IBBA behaviour. Box [k 1 ] is proved not to contain the minimizer k * because to Property (11), and similarly boxes [k 2 ] and [k 3 ] because to Property (12).
Through the creation of a subpaving, that is the union of non-overlapping boxes [30], the branch and bound method SIVIA provides an approximation of the feasible region of K. It bisects K in smaller boxes [k i ] until the constraint is proved to be fulfilled over [k i ] as a result of (14) or not to be fulfilled due to (11).
SIVIA algorithm stops when boxes [k i ] reach a minimum size . From Fig. 2

Proposal: GO as a SM tuning tool
This section is devoted to solve the problems of analysis and synthesis of SM using a GO approach. These problems are not convex in the general case. Hence, the global optimization algorithms presented just before are suited to solve such problems. Although the complexity of these algorithms grows exponentially with the number of variables, since they rely on a Branch and Bound structure, they are able to solve sliding mode problems of reasonable size as shown in the study-case. The main advantage of using global optimization over other optimization techniques (convex relaxation, local optimization, etc) to solve non-convex problems lies on the fact that the convergence to the global optimum is guaranteed.

SM analysis problem
Consider an objective sliding surface σ (x) = 0 satisfying the SM necessary condition L g σ = 0. The SM analysis problem involves verifying if u eq fulfills the existence condition (7). This problem can be posed as follows, sup δ∈Δ a(θ , δ) (15) with θ being a vector of constant tuning parameters given by the operator (for example coefficients of the sliding surface), δ the vector of variable parameters, Δ a subset of IR n δ with n δ the dimension of δ and a is the analysis function: IBBA can provide an enclosure [a, a] of the minimum a * useful to evaluate the SM condition over σ (x) = 0. If δ = x, then the equivalent control of the unperturbed system u * eq will be considered in (16) and the SM existence of the unperturbed system on a bounded region Δ of the state space will be evaluated in (15). That is, the nominal case will be addressed. On the contrary, if δ includes also all system uncertainties and external disturbances, then the equivalent control of the perturbed system u eq will be considered in (16) and the SM existence of the perturbed system will be evaluated in (15). That is, the robust problem will be addressed. From Property (10), we can derive Properties (17) and (18).
Property (17) provides a sufficient condition to prove that the system will slide on the sliding surface S over the subset Δ. Conversely, Property (18) provides a sufficient condition that the system will not slide over S in all Δ. In fact, the system will leave S at least at δ * the solution to Problem (15). When 0 ∈ [a, a], it is not possible to prove whether or not θ is a feasible solution since none of the conditions of Properties (17) and (18) is satisfied.

SM synthesis problem
The synthesis problem consists in choosing a sliding surface σ (x, θ) = 0 that guarantees the fulfillment of the SM existence condition for all admissible variations of system parameters and environment perturbations. This goal can be achieved by obtaining the set of all tuning parameters for which the SM existence condition holds, from where the designer could find θ by minimizing a cost function within this set. SIVIA algorithm and IBBA are suited to perform such computation.
Let Θ be a subset of IR n θ , l : R n θ → R be a cost function given by the system designer. It is assumed that an inclusion function of l is available, i.e. l has an analytic expression.
We define the analysis function at θ by: and is the minimum of Problem (15) with θ fixed. We also define the function The optimization problem (22) represents generically the synthesis task inf θ∈Θ l(θ ) The constraint of Problem (22) guarantees the SM existence condition and implies the resolution of the analysis problem (15). Using interval analysis, an enclosure of a sup over a box [θ ] can be provided [31]. Thus, IBBA can be used to solve Problem (22), while SIVIA can be used to characterize the set defined by the constraint In the literature, this kind of constraint is called a Semi Infinite Constraint (SIC), since it is equivalent to the infinite set of constraints Optimization problems involving SIC are called Semi Infinite Programs (SIP) and can be solved in a global way with different methods [32,33]. Furthermore, the characterization of the set defined by SICs has been studied in several works [34,35].

Application to the AUV heave direction control
The proposed technique is applied to the control of the AUV Ciscrea heave axis. Given the disturbances and uncertainties present in this type of application, the use of robust controllers is justified [36,37]. This problem will be analyzed from the controller design point of view, and series of tuning/analysis maps called subpavings will be generated.

AUV Ciscrea model
The model of the Ciscrea robot presented in Rosendo et al. [38], reduced to the heave direction, is used. This model is based on the ideas presented in [39] using parameters borrowed from [17]. Two coordinate systems are introduced: a NED-frame (North East Down) and a B-frame (Body fixed reference) for the localization as it is described in Fig. 1. The modeling of the heave direction is done by: where all the parameters and functions of the model are listed in Table 1  For more details about the mathematical model and numerical values used, refer to Rosendo et al. [38] where this model has been developed and experimentally tested. The AUV model (25)-(26) with the above considerations taken into account fits the nonlinear dynamical model (1) and the uncertainty in the control vector field is δg = −g · (ΔM R B + ΔM A ).

Global optimization for SM design
First, a desired SM dynamics is established. In particular, a first-order linear time response of the heave position is prescribed in which the time constant will be the tuning parameter to optimize. The sliding surface (3) is therefore defined as where z d is the reference position and λ is the tuning parameter. It can be easily proved that the necessary SM condition L g σ = 0 is fulfilled. The second step is to choose the Total propeller forces in heave direction minimization function. For the prescribed first-order dynamics, the optimal λ is chosen as the fastest one. Accordingly, the optimization function is l(θ) = −λ. The third step is to determine the expression of the equivalent control that is required by the optimization algorithm. The equivalent control for the unperturbed u * eq and perturbed system u eq are then derived: The limits u − and u + of the SM existence condition (7) are given by the maximum propeller torque |u − | = |u + | = τ max . For the AUV Ciscrea, τ max = 6 Nm in the nominal case. The discontinuous control is therefore τ pro = τ max sign(σ ).
Then, the optimization function, the available control effort, the expression of the equivalent control and the intervals of all variable parameters (including state variables) are fed to the optimization algorithm, which will return the optimal λ. This way, the SM of the system will be robust over the prescribed excursions of the state space and all other parameters and external disturbances.
Below, the variables selection for the optimization problem (22) is given: For the sake of brevity, the nominal model of the AUV was considered and just the environment torque τ env is considered as system disturbance. The ranges for speed and environment torque are consistent with the Ciscrea [38]. Particularly, note that the AUV speed is bounded by [−0.15, 0.15] m/s and the environment torque is limited to [−3, 3] Nm. Note also that no uncertainty in the model parameters is considered in this step, so just robustness against bounded environment torque is guaranteed. Of course, all uncertain parameters of the AUV model can be incorporated in δ, leading to a design robust against all these uncertainties too.
For the optimization problem posed above, the IBBA algorithm provides an enclosure [−0.3885, −0.3842] of the minimum. This means that λ = 0.3842 is the best feasible point, i.e. is the fastest λ that guarantees robust existence of SM over the variable intervals. If precision is increased, a bit faster λ ∈ [0.3842, 0.3885] could be obtained.

Comparison
It is interesting to compare the previous result with a traditional optimal tuning of SM [8]. The quadratic optimization technique applied to the AUV Ciscrea provides the optimal λ that minimizes the weighted quadratic errors in speed and position. For instance, -if position and speed are weighted equally, then the optimum is λ = 0.5 -if position is ten times more weighted than speed, then the optimum is λ = 0.31.
This optimization depends on the weight assignment, requires the system linearization and is completely decoupled from the SM existence conditions. In fact, after the sliding surface is designed, the control limits must be selected to guarantee SM for given excursions of variables and parameters or, alternatively, the admissible excursions are determined based on the available control effort. In contrast, the proposed optimization algorithm finds the optimal SM tuning subject to the available control effort and the admissible excursions of variables and parameters. A simulation example is presented in Fig. 3 where the results of a traditional SM manifold optimization and the proposed SM optimization subject to admissible variable excursions are compared. In Fig. 3a, the SM manifolds obtained from both optimizations are plotted with dotted lines, the red one corresponds to λ = 0.5 and the black one to λ = 0.38. The system was simulated with a constant disturbance τ tenv = −3, a maximum torque τ max = 6 and an initial position error e(0) = 0.5. The state trajectories are plotted with dashed lines. In black dashed line is shown the trajectory with the proposed design λ = 0.38. It is observed that the state starts from A and reaches the SM manifold at B. From then on, the state evolves along the prescribed manifold toward the equilibrium point. On the other hand, the state evolution corresponding to the standard optimization (λ = 0.5, red dashed line), reaches the optimal sliding manifold at point C, but the available control power is not enough to force the sliding regime. So, an undesirable openloop transient occurs until the surface is reached again at point D where the SM finally establishes. Therefore, for the standard optimal design to make sense, either the admissible state excursion should be reduced or the control power should be increased. In fact, e|(0)| < 0.2 is necessary to establish SM once the manifold is reached with τ max = 6 or, alternatively, τ max > 8.3 is required to establish SM once the manifold is reached from e(0) = 0.5.
The state trajectories for a smaller initial condition e(0) = 0.2 and for a higher control power τ max = 9 are shown with blue and green lines, respectively. Figure 3b shows the discontinuous control actions for all the cases discussed above using the same color code.

Global optimization for SM analysis and design assistance
The proposed SM design based on (22) finds the best tuning for the given optimization function, in this case study the fastest dynamics, subject to predetermined excursions of variables and parameters. Additionally, this work provides other powerful design and analysis tools. By means of (23), a series of maps can be developed providing much more insight into the robustness of the design. In fact, the generated maps help determining admissible excursions and operational limits, exploring the sensitivity of the optimal tuning with respect to state and parameter excursions, and manually tuning the controller. These maps show admissible excursions of variables and system parameters as function of design parameters. In the following, different optimization problems will be posed to illustrate the potentiality of the proposed tool. For the sake of clarity of presentation, 2D maps are presented. First, the subpaving λ versusż with the environment torque as variable parameter is built. This implies solving the problem given by (23), for the following variables and parameters: Figure 4 shows the resulting subpavings for three different torque magnitudes produced by the AUV thrusters, where the nominal case is at the center. These results are obtained applying SIVIA with = 0.01. In the subpavings, green boxes imply satisfaction of the imposed conditions, red boxes no satisfaction of them and blue boxes indicate that the algorithm cannot determine the fulfillment conditions. This methodology allows identifying in a guaranteed way where the objective dynamics will be achieved despite excursions in the parameters/states of the system. From the results, it is observed that the range of admissible dynamics increases as speed decreases (note area around theż = 0 axis). This behavior is due to the fact that at higher speeds the inertia of the AUV will require greater torque in order to guarantee the SM. Also it is noticed in Fig. 4 the existence of two lateral "branches" in the feasible zone, this is due to the non-linear behavior of the system (D N L ).

Remark 1
It is important to note that although the red boxes imply the non-satisfaction of the imposed conditions, this only means that there is at least a combination of the system parameters which SM is not guaranteed for.
It can be appreciated that a 25% reduction in torque (Fig. 4a) reduces the feasible area with respect to the nominal case (Fig. 4b), while an 25% increase widens it (Fig. 4a). It is possible to affirm from the results that even with variations of 25% in the available torque, the AUV will be able to work in the region of speedsż ∈ [−0.1, 0.1] with dynamics included in the interval λ ∈ (0, 0.38]. Figures 5 and 6 show the time response of the system and controller, while Fig. 7 depicts the state trajectories in the   Fig. 4b) and a larger one λ = 0.60 (outside the guaranteed area of Fig. 4b).
The systems starts from an initial position z = 2 while the set-point is z d = 1.5 An external disturbance τ env = −3 is applied from the beginning affecting the total torque. As expected, it is observed in Figs. 6 and 7 that the SM establishes immediately after the state reaches the sliding surface for λ = 0.38 and λ = 0.20. Figure 5 clearly shows that the time response with λ = 0.38 is much faster than with λ = 0.20. This means that choosing λ = 0.20 implies Fig. 7 State evolution in the phase space over-sizing the thrusters, or alternatively that the system is made unnecessarily slow. On the other hand, when λ = 0.60, the SM is not established but after some undesirable long transient. This is because the SM existence condition is not fulfilled at the first approach to the surface. So, although the sliding regime is much faster, the overall response in Fig. 5 is dominated by the open-loop reaching dynamics and even an overshoot is observed.
So far, the behavior of the system under external disturbances and variations in the control power has been analyzed. In the same way, it is possible to analyze what the effects of parametric uncertainties are. As the first case of analysis, Fig.  8 shows a variation of 5%, 10% and 25% in the M a parameter, observing a narrowing of the area with working conditions. For this case, referring to (23), the variables and parameters are:  (31) where the interval for M A corresponds to the ±5% excursion.
Obviously, the proposed methodology can be used to evaluate the effects of combined system uncertainties. Figure 9 shows the effect of a 25% simultaneous variation in M a and D N L parameters. It is possible to observe that the effect of D N L variation is dominant, considerably reducing the lateral branches (in comparison with Fig. 8). The optimizations problem in this case is posed as Another problem of interest for a designer is to determine the minimum control effort necessary to guarantee a given dynamics (λ value). This involves performing the subpaving of λ and |τ max | and solving (23) with: The resulting subpaving is shown in Fig. 10, which allows concluding that the fastest dynamics for τ max = 6 is λ = 0.38. As expected, this result coincides with the result of the first synthesis problem and Fig. 4b. Further, given a disturbance bound and a speed excursion, the fastest achievable dynamics (i.e. the largest λ) can also be identified. In other words, this allows determining the minimum control effort   (τ max ) needed to establish a prescribed convergence time for given intervals of state variables and uncertainties.
Finally, the procedure to determine the maximum disturbance that the system can endure without risking the SM existence is illustrated. With this aim, the subpaving with λ and τ env is performed. In this case (23) is solved for: Figure Fig. 11 shows the results for different control power. For the nominal case (Fig. 11b), it is concluded that the system works as designed within the area of interest even under a disturbance of |τ env | = 3. The subpavings for lower and higher control torque (Fig. 11a, c) show how sensitive the admissible external disturbance is to the control power.

Conclusions
A new technique for SM control tuning was devolved through the application of global optimization tools. The proposed approach returns a tuning that is optimal for a given criterion subjected to admissible excursions of state variables. Moreover, a robust design is achieved by incorporating also intervals for uncertain parameters. The joint consideration of SM dynamics, SM existence and uncertainties in the optimization procedure is the main distinctive feature of the proposal. Other advantages are that the optimization is global and does not require system linearization as other optimization techniques. Furthermore, the proposal provides additional tools for SM analysis and design by the user. In fact, by means of a subpaving analysis, the operational region, the sensitivity of the design to different uncertainties and the manual selection of tuning gains can be evaluated.
The potentiality of the proposal is illustrated by means of an AUV control design as case study. The approach allowed finding the fastest SM dynamics for given excursions of the AUV speed and the environment torque. Furthermore, admissible parameter inaccuracies, perturbations and state excursions without risking the desired SM operation were obtained.
Future research will be devoted to customize the optimization functions including both manifold and control parameters so that the SM dynamics and the control power are simultaneously designed.