Modeling, Analysis, and Control of Mechanical Systems under Power Constraints

Mo deling, Analysis, and Con trol of Mec hanical Systems under P o w er Constrain ts Gork em Secer Contr ol Systems Design Dept., ROKETSAN Missile Industries Inc., A nkar a, T urkey. Abstract Significan t improv ements ha ve been ac hieved in motion control systems with the av ailabilit y of high speed p o wer switc hes and micro computers on the market. Even though motor driv ers are able to provide high torque con trol bandwidth under nominal conditions, they suffer from v arious physical constrain ts whic h degrade b oth output amplitude and bandwidth of torque control lo op. In this con text, peak p ow er limit of a pow er source, as one of those constrain ts, has not b een fully explored from the con trol persp ective so far. A conv entional and practical w ay of considering p eak p ow er limit in control systems is to mo del it as a trivial torque saturation derived from the allo wable torque at maximum speed satisfying the constraint. Ho wev er, this mo del is ov erly conserv ative leading to p oor closed lo op p erformance when actuators operate below their maximum sp eed. In this pap er, nov el wa ys of incorp orating p eak p ow er limits in to b oth classical and optimal controllers are presen ted upon a theoretical analysis rev ealing its effects on stabilit y and performance. Key wor ds: Motion Con trol, Dynamical Systems, Stabilit y Analysis of Nonlinear Systems, Describing F unction Analysis, Linear Matrix Inequalities, Constrained Optimal Con trol, Con trol Ly apuno v F unctions. 1 In tro duction Since op erational lifetime of mobile systems is limited b y capacit y of their onboard p ow er supply and b y en- ergy expenditure of actuators, running time can only b e maximized by combined impro vemen ts in b oth supply side and actuation. In this con text, for the p ow er sup- ply , while some systems use energy harvesting mec ha- nisms [1–3], there is an interesting work [4] whic h adds autonomous rec harging capabilities to a mobile rob ot for prolonged service. F rom the viewp oint of actuation, energy consumption can be minimized either by seeking optimal mechanical design [5–7] for a particular behav- ior or by designing energy efficient controls [8, 9]. F o cus of this pap er is on the latter, particularly on motion con- trollers resp ecting instantaneous p ow er supply limits. Motion con trol systems consist of actuators and driv ers. While the actuator generates torques/forces necessary to mo v e the system and/or payload attached to it, the driv er is resp onsible for controlling the actuator. In this con text, electrical motor-based systems are b ecoming extremely p opular thanks to impro v emen ts in mag- net materials (e.g., NdF eB) and in v erter tec hnology . P ermanen t-magnet sync hronous motors (PMSMs) are highly preferred within this trend compared to other A C and DC machines due to their efficiency and torque densit y [10]. Drivers control PMSMs b y con v erting bus v oltages/curren ts to motor voltages/curren ts with ra- tios adjustable via dut y cycles of switc hing pulse width mo dulation signals. Developmen ts in high sp eed dig- ital devices and switc hing circuits hav e increased the p erformance of torque controllers on driv ers suc h that relativ ely high bandwidth can be pro vided allowing mo- tion control designers to ignore electrical dynamics and inefficiencies. Based on this abstraction, it is often suf- ficien t to model the driver and motor combination as a p erfect torque source. How ever, motion control systems suffer from v arious nonlinearities due to physical limi- tations in driving circuit and p ow er supply . There are t w o main nonlinearities present in motion controllers in this context: I) P eak curren t rating of switc hes (e.g., MOSFETs) whic h can b e modeled as saturation of motor torques u ∈ R n in the form sat( u, u max ) = [sat( u 1 ) sat( u 2 ) · · · sat( u n )] T (1) with sat( u i ) = sign( u i ) min {| u i | , u max } , i = 1 , 2 , . . . n . Due to saturation, actual and desired controller output b ecome different, th us leading controller states to b e up- dated wrongly . This may result in large ov ersho ots at plan t output and sometimes ev en instabilit y (see [11, 12]) when the controller has slo w and/or unstable mo des (i.e., Preprin t submitted to Journal of L A T E X T emplates No v em b er 15, 2021 p oles and ze ros close to and/or on the right hand side of imaginary axis, resp ectively). Numerous tec hniques and extensions to controllers, such as conditional integra- tion [13, 14], reference go vernors [15, 16], v ariable struc- ture con trollers [17] or anti windup with model recov ery filter [18], hav e been prop osed in the literature to handle saturation effects in terms of both stabilit y and p erfor- mance. Most of these tec hniques ha v e pro v able stabilit y results whic h are generally obtained b y describing func- tion analysis of saturation nonlinearit y [19], linear ma- trix inequalities [20, 21] or con vex programming [22, 23]. I I) Po wer limit of the onboard p ow er supply batteries and circuits [24], whic h is our focus in this pap er. Even though this constraint takes place in electrical part of the system, it can b e mapped to motion control level b y enco ding the constraint into softw are as psat( u, ˙ q ) = [psat( u 1 , ˙ q 1 , ¯ P 1 ) · · · psat( u n , ˙ q n , ¯ P n )] T where psat( u i , ˙ q i , ¯ P i ) =      u i if P i ≤ ¯ P i − ˙ q i ± q ˙ q 2 i + 4 ¯ P i R i / ( k t ) 2 i 2 R i / ( k t ) 2 i else if u i ≷ 0 , ˙ q i is the generalized velocity of the i th join t, ¯ P i is the p ow er budget allo cated to i th join t supplied b y a common source ha ving peak pow er P max with P max = P n i =1 ¯ P i , and the p ow er input P i to the i th motor includes mec hanical pow er output determined b y torque constan t ( k t ) i and motor losses dissipated across the resistance R i as P i := u i ˙ q i + u 2 i R i / ( k t ) 2 i . (2) F or high efficiency motors, electrical losses can b e ig- nored whic h leads to a simpler expression for p ow er limit saturation as psat( u i , ˙ q i ) =  u i if P i ≤ ¯ P i ¯ P i / ˙ q i otherwise. (3) The p o w er limit is a nonlinearity , which might lead to instabilit y and deterioration in closed loop p erformance. In con trary to the extensive amount of w ork on torque saturation, we ha ve not found any w ork considering the external p ow er supply limit constraint explicitly to this date, despite the fact that it has become particularly im- p ortan t for actuators with recen t adv ances in mobile sys- tems. A relativ ely reasonable yet conserv ative appro xi- mation that might b e used in practice is psat( u i , ˙ q i , ¯ P i ) ≈ sat( u i , ¯ P i ¯ v i ) (4) where ¯ v i denotes the no-load v elo city of a motor located at i th join t. How ev er, this may result in sup erfluous ac- tuator and pow er supply selections to meet p erformance sp ecifications of the motion controller, even tually limit- ing the ov erall system p erformance esp ecially in the ap- plications such as legged rob ots or electric v ehicles where total w eight is crucial for energetic efficiency . T o sho w that, exact mo del of the p ow er supply limit nonlinear- it y is analyzed in this paper from v arious persp ectives in comparison with the approximate mo del (4). In particu- lar, we deriv e the describing function of the p ow er limit nonlinearit y and maxim um closed lo op frequency attain- able in its presence to understand gain and phase effects on stability and steady-state p erformance, resp ectively . This analysis is then used to motiv ate design of stable con trollers based on exact mo del of the nonlinearity for Euler-Lagrange mec hanical systems. These controllers require that the total p ow er limit of the battery is allo- cated to individual p ow er budget of each actuator unit (i.e., ¯ P i ) for whic h a static solution strategy (i.e., assign- ing a fixed p o w er limit prior to operation) is used con- v en tionally . Note that, the domain of p ow er saturation function can b e reduced to psat( u, ˙ q ) = [psat( u 1 , ˙ q 1 ) · · · psat( u n , ˙ q n )] T . when the static allo cation is employ ed. On the other hand, p erformance can b e improv ed by up dating p ow er budget dynamically with time-v arying limits ¯ P i ( t, q i , ˙ q i ). F or example, when some join ts consume less pow er than their limits, c hannelling their lefto ver p ow er to remain- ing join ts migh t improv e the p erformance. T o maximize this adv antage, a general control arc hitecture com bining dynamic p o w er allo cation with exact p o w er limit model is developed for linear Euler-Lagrange mechanical sys- tems in this paper. The prop osed architecture is p osed as an optimization problem to solve a finite horizon op- timal control problem and dynamic p ow er allo cation si- m ultaneously . Note that, unless otherwise sp ecified, the static allo cation scheme is used throughout the pap er. The organization of this paper is as follows : Section 2 presen ts frequency domain analysis of the p ow er limit nonlinearit y that pro vides motiv ation to design con- trollers accoun ting for this nonlinearity . After such classical con trollers are proposed in Section 3, they are extended to a new class of algorithms in Section 4 such that dynamic p ow er allocation and feedbac k control is unified under a general architecture. Section 5 gives preliminary exp erimen tal results on a single degree-of- freedom (DOF) actuator prototype. Finally , conclusions and future work are presented in Sec. 6. 2 F requency Response Characteristics In this section, effects of pow er limit psat (3) (i.e., the simpler model without electrical losses) on closed-lo op 2 stabilit y and p erformance are inv estigated in compar- ison to the approximate mo del (4). In particular, the sin usoidal input describing f unction of nonlinearities (i.e., a first harmonic approximation to the nonlinearity) are derived to study destabilizing c haracteristic of b oth phenomenons, whereas the maximum closed loop band- width is derived to determine a b enchmark for linear con trollers using different strategies for p ow er supply limit. Findings will b e used as foundational premises to motiv ate more in-depth exploration of using exact mo dels of p o w er supply limit for con troller design in the subsequen t section. 2.1 Describing F unction A nalysis Describing function of a nonlinearit y is a frequency do- main approximation known as quasi-linerization [25]. Ev en though it was introduced long ago, it is still widely used to predict oscillations of linear dynamical systems sub jected to a nonlinearit y [26–28]. In this con text, the nonlinearit y is approximated b y its fundamental har- monic frequency resp onse under a sinusoidal input to iden tify limit cycles and study their local behavior. This appro ximation often pro vides sufficient accuracy for as- sessmen t of p erio dic b eha viors, mo deled b y frequency domain resp onse, in terms of stability and gain/phase c haracteristics since linear systems attenuates higher or- der harmonics considerably due to their low-pass filter c haracteristics [29]. In particular, the main fo cus of this pap er is on motion con trol of electromec hanical systems whose dynamics hav e at least 2 nd order relative degree, hence filtering out higher-order harmonics to a large ex- ten t. Consider a single DoF actuator whose transfer function from torque u to velocity ˙ q is denoted b y G ( s ). Denot- ing the describing function of the p ow er saturation b y N ( A, ω ), transfer function and describing function can b e expressed in phasor notation as G ( s ) = X G ( ω ) φ G ( ω ) N ( A, ω ) = X N ( A, ω ) φ N ( A, ω ) . Cascade connection of the describing function and trans- fer function is defined by G ( s ) N ( A, ω ) = X ( A, ω ) φ ( A, ω ) with X ( A, ω ) = X G ( ω ) X N ( A, ω ) φ ( A, ω ) = φ G ( ω ) + φ N ( A, ω ) . First, consider the case where p ow er saturation is not active, i.e., P ≤ P max as defined in (3). Since u = psat( u, ˙ q ), the describing function b ecomes N ( A, ω ) = 1. F or the other case with P > P max , we start the deriv ation of describing function by assuming a sinusoidal torque input u = A sin( ω t ) In this case, the v elo cit y and required p ow er take the form ˙ q = AX sin( ψ + φ ) P = A 2 X sin( ψ + φ ) sin ψ = A 2 X 2 [cos φ − cos(2 ψ + φ )] with ψ := ω t . Observe that the required p ow er P is twice the angular frequency ω of the excitation torque, mean- ing that it is sufficien t to consider the scaled time quan- tit y in the in terv al ψ ∈ [0 , π ] for the following analysis. Since required p ow er is a contin uous function of ψ , in- v erse image of the open interv al Ψ = { ψ : P ( ψ ) > P max } is op en [30], i.e., P − 1 (Ψ) ∼ = ( ψ l , ψ u ) ⊂ [0 , π ] with P ( ψ l ) = P ( ψ u ) = P max . (5) Solutions to low er and upp er bounds in (5) are obtained as ψ ( l u ) = ± arccos(cos φ − 2 P max / ( A 2 X )) − φ 2 mo d π . Therefore, when ψ ∈ ( ψ l , ψ u ) (i.e., active region of p o w er supply limit), pow er saturation maps desired motor torques to ¯ u := psat( u, ˙ q ) = P max / ( AX sin( ψ + φ )) . No w that input-output torque relations of pow er supply limit are obtained in b oth cases, F ourier co efficients of the p ow er limit nonlinearity o ver a p erio d of the torque input u = A sin( ω t ) i.e., t ∈ [0 , 2 π ], can b e computed as c N := X N cos φ N = 1 π A Z 2 π 0 psat( u, ˙ q ) sin ψ d ψ s N := X N sin φ N = 1 π A Z 2 π 0 psat( u, ˙ q ) cos ψ d ψ. (6) Defining the follo wing notation for limits of the integral, Z u 1 ··· u n l 1 ··· l n h ( x )d x = n X i =1 Z u i l i h ( x )d x, in tegrals in (6) can b e partitioned into active and passiv e regions of p ow er saturation, corresponding to P ≤ P max and P > P max , resp ectively , as c N = 2 π A    ψ l ,π Z 0 ,ψ u A sin 2 ψ d ψ + ψ u Z ψ l ¯ u sin ψ d ψ    s N = 1 π A    ψ l ,π + ψ l , 2 π Z 0 ,ψ u ,π + ψ u A 2 sin 2 ψ d ψ + ψ u ,π + ψ u Z ψ l ,π + ψ l ¯ u cos ψ d ψ    . 3 After some calculus to ev aluate integrals, the following nonlinear system of equations are obtained c N = 2 P max π A 2 X Y N + 2( π − ∆ ψ ) + sin 2 ψ u − sin 2 ψ l 2 π s N = 2 P max π A 2 X Z N + cos 2 ψ u − cos 2 ψ l 2 π with L ψ := ln sin( ψ u + φ ) sin( ψ l + φ ) , Y N := ∆ ψ cos φ − L ψ sin φ , Z N := ∆ ψ sin φ + L ψ cos φ , and ∆ ψ := ψ u − ψ l . Since these equations do not admit analytical closed form solutions, n umerical metho ds can b e employ ed to obtain solution to X N and φ N at a given amplitude A and frequency ω , which define the describing function N ( A, ω ) = X N ( A, ω ) φ N ( A, ω ) . Example 1 Consider a p osition servo system char ac- terize d by a se c ond-or der close d-lo op tr ansfer function T ( s ) = ω 2 n s 2 + 2 ζ ω n s + ω 2 n with natur al fr e quency ω n = 50 π r ad/s and damping r atio ζ = 0 . 8 . Me chanics of the servo system c onsists of a mass m = 1 k g m 2 and visc ous damping d = 0 . 05 N s/r ad c orr esp onding to a plant with tr ansfer function G ( s ) = 1 sm + d . (7) The system is c ontr ol le d by a pr op ortional-derivative (PD) c ontr ol ler in the form u = K p ( q d − q ) − K d ˙ q with desir e d p osition q d and PD gains K p = ω 2 n and K d = 2 ζ ω n − d , r esp e ctively. Assume ther e exists a p ower limit P max = 400 W . We c onsider b oth exact and appr oximate mo dels of the p ower limit to evaluate their effe cts on the c ontr ol lo op thr ough describing functions c ompute d for amplitudes in the r ange A ∈ [1 , 500] r ad . Note that the appr oximate mo del is establishe d as a standar d tor que satur ation u max = 100 N m which fol lows fr om the assumption that sp e e d is limite d by ˙ q max = 4 r ad/s . T o b e gin with, the op en-lo op tr ansfer function of the PD c ontr ol le d me chanic al system including a describing func- tion is evaluate d on a Nyquist plot in Fig. 1. It is ob- serve d that magnitudes of b oth nonline arity mo dels ar e less than or e qual to unity, and the exact mo del pr o duc es non-ne gative imaginary p art wher e as appr oximation r e- sults in pur ely r e al describing function. These findings le ad to two main c onclusions: 1 ) Nonline arities have no dir e ct effe ct on stability sinc e their describing functions do neither interse ct with or encir cle d by the op en-lo op curve. 2) Nonline arities incr e ase the phase mar gin c om- monly, henc e the r obustness of the close d lo op system. However, it should b e distinguishe d that exact nonline ar- ity mo del pr ovides mor e phase mar gin (even mor e with incr e ase in amplitude) sinc e it brings phase le ad into the system b e c ause of its p ositive imaginary p art as opp ose d to satur ation. Figure 1. Nyquist diagram of nonlinearity describing func- tion (dotted green), op en-lo op transfer function in the ab- sence of pow er constraint (solid red), and open-lo op transfer function when coupled with nonlinearity (dashed blue). 2.2 Maximum Close d L o op Bandwidth In order to further impro ve on results of the last subsec- tion summarized in Figure 1, the closed loop bandwidth is ev aluated more rigorously in this subsection by con- sidering other nonlinearities common in motion control systems. While doing so, the exact mo del of pow er limit is quantitativ ely compared to its approximate mo del in terms of the achiev able maximum bandwidth in the lin- ear region of dynamics. Therefore, describing function of nonlinearity mo dels are not needed. Consider a p osition con troller with a sin usoidal reference q d = Y cos( ω t − θ ) applied to a single DoF linear system m ¨ q = u − k q − d ˙ q − τ c sign( ˙ q ) with inertia/mass m , viscous friction coefficient d , and Coulom b friction τ c . W e are interested in maximum attainable bandwidth of the closed-lo op system in the linear region of dynamics under different ph ysical con- strain ts. First, consider a p ow er supply limit P max (i.e., P ≤ P max ) and no-load speed limit − ˙ q max ≤ ˙ q ≤ ˙ q max . F or the remaining cases, standard torque saturation − u max ≤ u ≤ u max with the conserv ative approxima- tion to torque limit, i.e., u max = P max / ˙ q max giv en in (4), is mo deled in addition to the no-load sp eed limit. A t cut-off frequency ω c , tracking the reference signal as q = ¯ Y cos( ω c t ) with am plitude ¯ Y := Y / √ 2 and phase dela y θ requires velocity , torque and p ow er ˙ q = − ¯ Y ω c sin( ω c t ) u = ¯ Y  ( k − mω 2 c ) cos( ω c t ) − d ω c sin( ω c t )  − τ c sign(sin( ω c t )) P = − ¯ Y 2 ω c 2  ( k − mω 2 c ) sin(2 ω c t ) + d ω c (cos(2 ω c t ) − 1)  − ¯ Y ω c τ c | sin( ω c t ) | . In order to c haracterize the maxim um closed loop band- width under differen t constrain ts, a common n umerical 4 optimization problem is formulated as follows : maximize ω c ω c sub ject to ω c > 0 · · · (1) max 0 ≤ t ≤ 2 π/ω c | ˙ q ( t ) | ≤ ˙ q max · · · (2) max 0 ≤ t ≤ π/ω c P ( t ) ≤ P max · · · (3) max 0 ≤ t ≤ 2 π/ω c | u ( t ) | ≤ u max · · · (4) F or the first case, whic h uses the exact mo del of pow er supply limit, inequality (4) is remo ved whereas inequal- it y (3) is remo v ed for the remaining case (i.e., appro xi- mate mo del). Example 2 Consider the same setup use d in the Nu- meric al Example of Se c. 2.1 (i.e., r otational me chanic al system, motor, and the driver). The the or etic al p ermis- sible b andwidth of the system with exact and appr oxi- mate mo dels ar e evaluate d for amplitudes in the r ange Y ∈ (0 , 1] . This analysis is r ep e ate d for differ ent p ow er supply limits, sp e cific al ly P max ∈ { 200 , 400 , 600 } W . R e- sults ar e il lustr ate d in Fig. 2 by a r atio of p ermissible b andwidths of the exact mo del to the appr oximate mo del, denote d by ω psat c and ω sat c , r esp e ctively. Using exact mo del incr e ases the maximum close d lo op fr e quency substan- tial ly, ne arly doubling the b andwidth at low amplitude r e- gion (e.g., Y < 10 de gr e es for P max = 600 W ). F urther- mor e, it is se en fr om the figur e that the r atio, interpr ete d as a function of amplitude Y , is not differ entiable at two p oints. This c an b e explaine d by the change of active non- line arity. In p articular, the sp e e d limit, which is p assive for b oth mo dels at smal l amplitudes, b e c omes active first for the exact mo del with the inactivation of p ower limit, and then for the appr oximate mo del with the inactiva- tion of tor que limit. F or higher amplitudes, b oth systems ar e c onstr aine d by the sp e e d limit alone, which makes the b andwidth r atio unity. 0 0.2 0.4 0.6 0.8 1 0.5 1 1.5 2 2.5 Figure 2. Maxim um closed-lo op frequency as a function of input amplitude for different springs. 3 Con troller Design with Static P ow er Allo ca- tion F requency domain analysis in Sec. 2 assesses p erfor- mance and lo cal stability of an actuator supplied b y an energy source with a pow er supply limit in comparison to an actuator under a torque saturation corresp ond- ing to the conserv ative approximation (4) of the pow er limit. In particular, p erformance is measured b y theo- retical bandwidth limit of an actuator in the linear op- eration region, whereas describing function can serve as a to ol to conduct lo cal stability analysis around pre- dicted limit cycles. In summary , the presen ted analysis suggests that using an exact mo del of the p ow er supply limit outperforms the appro ximate mo del in terms of the op en-lo op phase margin, whic h is an indicator of closed- lo op robustness, and closed-lo op bandwidth. These re- sults motiv ate further in vestigation of exact mo dels for the pow er supply limit, especially during con troller de- sign. This section explicitly fo cuses on the wa ys of in- cluding p ow er limit into existing motion con trol frame- w orks. In particular, w e consider classical con trollers for Euler-Lagrange mechanical systems. 3.1 Plant Linear and nonlinear mec hanical systems are studied. Consider first, linear systems. Underactuated and fully- actuated linear mechanical systems can be commonly describ ed by M ¨ q + D ˙ q + K q = S psat( u, ˙ q ) (8) with mass matrix M = M T > 0, friction D = D T ≥ 0, stiffness K = K T ≥ 0, and actua- tor selection matrix S . F ollowing, [31], the generalized co ordinates can b e decomp osed into passive (unactu- ated) and actuated join ts as q =  q T u q T a  T with unac- tuated join t vector q u ∈ R n u and actuated join t vector q a ∈ R n a , resp ectively . Hence, the actuator selection matrix takes the form S = " 0 n u × n a I n a × n a # . On the other hand, fully-actuated rigid manipulators with nonlinear dynamics take the form of Euler Lagrange systems with Ra yleigh dissipation [32], whic h admit dy- namics M ( q ) ¨ q + C ( q , ˙ q ) ˙ q + D ˙ q + G ( q ) = psat( u, ˙ q ) (9) with mass matrix M ( q ) = M T ( q ) > 0, Coriolis ma- trix C ( q , ˙ q ), gravitational force G ( q ) derived from p o- ten tial field U ( q ) for all generalized co ordinates and ve- lo cities in tangent bundle of configuration space Q , i.e., 5 ∀ ( q T , ˙ q T ) T ∈ T Q . Note that Coriolis matrix satisfies the sk ew-symmetric prop erty [33] ˙ q T h 1 2 ˙ M ( q ) − C ( q, ˙ q ) i ˙ q = 0 , ∀ q , ˙ q ∈ R n . (10) 3.2 Contr ols 3.2.1 Line ar Systems F or the plan t with linear dynamics in (8), assume a linear dynamic controller of the form ˙ x c = A c x c + B p q + B d ˙ q u = C x c + K p q + K d ˙ q (11) to stabilize the system in the absence of p ow er limit sat- uration where x c ∈ R n c . F urthermore, the controller can b e augmented b y adding anti-windup comp ensation to a v oid windup of con troller states x c for which we con- sider b oth static and dynamic approac hes. In particular, conditional integration (CI) and mo dern anti-windup (MA W) techniques [18] with a general form ulation are ev aluated. These techniques are selected b ecause of the follo wing reasons : 1) CI is commonly used in practice since it is simple in the sense that it do es not in tro- duce any additional parameters for anti-windup aug- men tation [34]. 2) MA W provides superior p erformance compared to static techniques b y guaranteeing stabilit y and optimal p erformance with the additional flexibil- it y enco ded by the dynamics of anti-windup augmenta- tion [35]. W e prop ose to up date the controller states according to the CI rule ˙ x c = N [ ρ ( C, H )] ( A c x c + B p q + B d ˙ q ) (12) where H = { 1 ≤ i ≤ n a : P i > ¯ P i } denotes the set of actuator channels demanding more p o w er than the limit, ρ ( C, H ) selects the rows of C b elonging to the set H ρ ( C, H ) = [ c i ] , ∀ i ∈ H with ρ ( C , ∅ ) = I , and N [ ρ ( C , H )] denotes the pro jection op erator onto n ull-space of the matrix ρ ( C , H ) formu- lated as N [ ρ ( C, H )] = V H Σ H ( V H ) T (13) with singular v alue decomp osition ρ ( C, H ) = U H Σ H V T H . The null-space pro jection freezes the integrator of state v ector in the subspace associated with saturating actua- tor c hannels. Under the controller action, the closed loop system can b e written as ˙ x = A ( H ) x + B σ ( x ) (14) where x :=     q ˙ q x c     , σ ( x ) := psat( u, ˙ q ) − u, u = κx, A ( H ) :=     0 I 0 M − 1 ( K p + K ) M − 1 ( K d + D ) M − 1 C N [ ρ ( C, H )] B p N [ ρ ( C, H )] B d N [ ρ ( C, H )] A c     , B :=     0 M − 1 S 0     , and κ := h K p K d C i . As a second anti-windup alternative, MA W scheme can b e used in which con troller dynamics and con troller out- put take the form ˙ x c = A c x c + B p q + B d ˙ q + E c σ ( x ) u = κx + E σ ( x ) (15) where x c no w includes an ti-windup augmen tation states, and E c and E denote an ti-windup compensation gains con tributing to the controller dynamics and output, re- sp ectiv ely . Thus, the closed lo op system dynamics in (14) is obtained with A ( H ) :=     0 I 0 M − 1 ( K p + K ) M − 1 ( K d + D ) M − 1 C B p B d A c     and B :=     0 M − 1 S ( I + E ) E c     . F or the function σ ( x ) whic h is common in b oth con troller structures, a closed conv ex p olytop e can b e defined as L = { x ∈ T Q × R n c : | γ i [ K c ] i x | ≤ ¯ P i / ¯ v i , i = 1 , . . . , n a } with [ K c ] i = ρ ( K c , { i } ) and a p ositive scalar γ i < 1, suc h that ∀ x ∈ L , the following holds 0 ≤ γ i [ K c ] i x ≤ psat([ K c ] i x, ˙ q ) ≤ [ K c ] i x if [ K c ] i x ≥ 0 [ K c ] i x ≤ psat([ K c ] i x, ˙ q ) ≤ γ i [ K c ] i x ≤ 0 otherwise . This translates to − (1 − γ i )[ K c ] i x ≤ [ σ ( x )] i ≤ 0 if [ K c ] i x ≥ 0 0 ≤ [ σ ( x )] i ≤ − (1 − γ i )[ K c ] i x otherwise . (16) 6 Theorem 3 The system in (14) is asymptotic al ly stable in L if x T h ( A ( H ) + B Π( H )) T Q + Q ( A ( H ) + B Π( H )) i x < 0 for al l H ∈ P ( n a ) with a matrix Q = Q T > 0 , P ( n a ) denoting the p ower set of inte gers fr om 1 to n a , which, having c ar dinality of 2 n a , δ ( H ) a mapping fr om indic es of satur ating joints to R n × n c define d as Π( H ) =        − δ (1 , H )(1 − γ 1 )[ K c ] 1 − δ (2 , H )(1 − γ 2 )[ K c ] 2 . . . − δ ( n a , H )(1 − γ n a )[ K c ] n a        , and a binary indic ator function δ ( i, H ) =  1 if i ∈ H 0 otherwise. Pro of 1 Consider a Lyapunov function V ( x ) = x T Qx . The set L is asymptotic al ly stable if ˙ V ( x ) = x T  A ( H ) T Q + QA ( H )  x + 2 x T QB σ ( x ) < 0 , (17) for al l x ∈ L . In other wor ds, it c an b e said that the set L is inside the r e gion of attr action (RO A). Sinc e ˙ V ( x ) is line ar in σ ( x ) and A , it attains its maximum in L at one of the vertic es define d in (16), which implies ˙ V ( x ) ≤ max H∈ P x T  A ( H ) T Q + QA ( H )  x +2 x T QB Π( H ) x. Henc e, (17) is satisfie d if, for al l H ∈ P , we have x T h ( A ( H ) + B Π( H )) T Q + Q ( A ( H ) + B Π( H )) i x < 0 . The RO A can be estimated b y fitting a prescribed shape to the interior of p olytop e L [22] and [36]. In this pa- p er, ellipsoid E ( Q, α ) = { x : x T Qx ≤ α } ⊂ L , identi- fied by the Lyapuno v function and α > 0, is chosen as the reference shap e. Therefore, volume of RO A can b e maximized b y join tly tuning the scalar γ , which deter- mines the affine inequalities of p olytop e L , and α , which determines the b oundary of ellipsoid, in line with prior w orks [23] and [36]. Using Sch ur’s complement, this can b e translated to a constrained optimization problem de- fined in the language of LMIs as max A c ,B p ,B d C,K p ,K d W, [ γ i ] log det W s.t. W > 0 , 0 < γ i < 1 ∀ i = 1 . . . n a , W ( A ( H ) + B Π( H )) T + ( A ( H ) + B Π( H )) W < 0 ∀H ∈ P ( H ) , γ 2 i h i W h T i ≤ 1 ⇔ " 1 γ i h i W γ i W h T i W # ≥ 0 ∀ i = 1 . . . n a (18) with W :=  Q α  − 1 and h i := ¯ v i ¯ P i [ K c ] i . Unfortunately , this program is not conv ex for join t optimization of con- troller gains and remaining parameters, whic h are W and γ i s, to find the controller with maximal RO A since they are not simultaneously linear in the LMIs. Itera- tiv e methods can b e used to decomp ose the problem in to con v ex subproblems and solv e each stage of nested op- timization via conv ex metho ds [36, 37]. Nevertheless, in Sec. 5, a controller with a CI an ti-windup augmen tation based on the exact mo del of p o w er supply limit is de- signed following the pro cedure developed so far in this section and ev aluated exp erimentally in b oth time and frequency domains compared to the same controller us- ing the approximate mo del (3) of the p ow er limit. 3.2.2 Nonline ar Systems Consider no w a fully-actuated nonlinear Euler-Lagrange system controlled b y a passivity-based static nonlinear con troller (also kno wn as PD plus gravit y compensation) in the form u = G ( q ) − K p q − K d ˙ q . (19) with K p = K T p > 0. A Lyapuno v function candidate V ( q, ˙ q ) = 1 2 ˙ q T M ( q ) ˙ q + 1 2 q T K p q , (20) can b e derived from the kinetic energy along with the ”shap ed” p otential energy [38] [39]. This Ly apunov func- tion candidate is radially unbounded and p ositive defi- nite. Its deriv ative along the tra jectories of (9) can b e found as follows : First differentiate it with resp ect to time ˙ V ( q, ˙ q ) = ˙ q T M ( q ) ¨ q + 1 2 ˙ q T ˙ M ( q ) ˙ q + ˙ q T K p q . (21) Then, substituting system dynamics in (9) gives ˙ V ( q, ˙ q ) = ˙ q T (psat( u, ˙ q ) − C ( q , ˙ q ) ˙ q − D ˙ q − G ( q )) + 1 2 ˙ q T ˙ M ( q ) ˙ q + ˙ q T K p q (22) 7 whic h can b e further simplified to ˙ V ( q, ˙ q ) = ˙ q T psat( u, ˙ q ) + 1 2 ˙ q T  ˙ M ( q ) − 2 C ( q, ˙ q )  ˙ q + ˙ q T K p q − ˙ q T G ( q ) − ˙ q T D ˙ q. (23) Using the skew-symmetric prop erty in (10) and substi- tuting the con trol action in (19) after adding/subtracting the term ˙ q T K d ˙ q , we obtain ˙ V ( q, ˙ q ) = ˙ q T psat( u, ˙ q ) − ˙ q T u − ˙ q T ( K d + D ) ˙ q . (24) Consider, first, P i ≤ ¯ P i for all join ts meaning that psat( u, ˙ q ) = u . This yields ˙ V ( q, ˙ q ) = − ˙ q T ( K d + D ) ˙ q ≤ 0 (25) negativ e as long as ˙ q 6 = 0 for all system tra jectories. With psat( u, ˙ q ) = u , the system dynamics reduce to M ( q ) ¨ q + C ( q , ˙ q ) ˙ q + K p q + ( K d + D ) ˙ q = 0 (26) whic h renders the origin ( q = 0 and ˙ q ) unique equi- librium p oint in the state space. Then, application of LaSalle’s In v ariance Principle [40] ensures that the ori- gin is globally asymptotically stable (GAS). F or the re- maining case with P i > ¯ P i , the time-deriv ative of the candidate Lyapuno v function becomes strictly negativ e for all q and ˙ q with ˙ V ( q, ˙ q ) = − ˙ q T ( K d + D ) ˙ q − X i ∈ F P i − ¯ P i < 0 , (27) th us concluding that the system under passivity-based con trol action sub jected to p ow er supply limits is GAS. T able 1 Ph ysical parameters of the tw o-link rob ot mo del. P arameter Sym b ol V alue Link 1 mass m 1 16 kg Link 2 mass m 2 12 kg Link 1 inertia I 1 18 kg m 2 Link 2 inertia I 2 7.5 kg m 2 Join t damping matrix D diag(10,10) Nm s / rad Link 1 length h 1 1 m Link 2 length h 2 1 m Link 1 angle q 1 ( − ) rad Link 2 angle q 2 ( − ) rad Link 1 torque u 1 ( − ) N m Link 2 torque u 2 ( − ) N m Gra vit y g 9.8 m / s 2 Example 4 Consider a manipulator il lustr ate d in Fig. 3. Assume that joint actuators ar e supplie d by a single p ower sour c e with p e ak p ower P max = 2 k W , c orr e- sp onding to p ower limit of e ach joint ¯ P 1 = ¯ P 2 = P max / 2 . Using p ar ameter values given in T able 1 and a p assivity- b ase d c ontr ol ler K p = diag  m 1 ω 2 n , m 2 ω 2 n  K d = diag ([2 m 1 ζ ω n , 2 m 2 ζ ω n ]) , with ω n = 2 π √ 2 and ζ = 0 . 9 , the r ob ot is simulate d fr om initial c onditions q 1 = − π/ 2 and q 2 = π . As il lustr ate d in Fig. 4, joint tr aje ctories and joint tor ques verify that the origin is asymptotic al ly stable for this system even though me chanic al p ower of e ach joint actuator is upp er b ounde d by the p ower supply limit of 1 kW . 4 A Unified Algorithm for Dynamic Allocation and Con trol Classical controllers with static p ow er allocation enforce P j ≤ ¯ P j for all joints. How ever, this may be unneces- sary in some cases b ecause channelling more p ow er to particular joints (i.e. P j > ¯ P j ) without compromising P n a i =1 P i < P max is sometimes a feasible strategy which can improv e tracking p erformance. F or example, when a join t requires pow er less than its predefined limit for a desired tra jectory , its residual pow er ma y be used tem- p orarily by another joint which uses its av ailable p ow er completely but still fails to follo w its o wn desired tra jec- tory . This motiv ates the use of dynamic p ow er allo ca- tion for improv ed p erformance. In the following subsec- tions, tw o quadratic programming (QP) problems are presen ted to sim ultaneously solve con trol and dynamic p o w er allo cation problems for linear and nonlinear sys- tems formulated as a finite-horizon optimal control prob- lem and a control Lyapuno v function-based QP prob- lem, resp ectively . 4.1 Line ar Systems Without loss of generalit y , we consider linear mechani- cal systems in (8) extended with constan t gra vitational Figure 3. F ully-actuated tw o-link planar manipulator. 8 -1 0 1 2 3 -2 -1 0 1 0 0.5 1 1.5 -4 -2 0 2 4 Figure 4. T ra jectories (top), torques (middle), and p ow er consumption (b ottom) of a 2-DoF rotary rob ot with a pas- sivit y-based con troller. disturbance G . This results in state-space mo del ˙ x = F c x + H c υ + g c (28) with F c = " 0 I − M − 1 K − M − 1 D # , H c = " 0 M − 1 S # , x = " q ˙ q # , and g c = " 0 − M − 1 G # . Zero-order-hold discretization of dynamics with sample time ∆ t yields x k +1 = F x k + H υ k + g (29) with H =   ∆ t Z 0 e F c t d t   H c , g =   ∆ t Z 0 e F c t d t   g c , and F = e F c ∆ t . Giv en an initial condition x 0 and a sequence of con trol action { υ i } , the solution to discrete dynamics is obtained as x n = F n x 0 + n − 1 X i =0 F n − 1 − i ( H υ i + g ) . (30) W e are interested in finite-horizon optimal control of dy- namical system (29) in the presence of linear constrain ts on state and con trol inputs and p ow er supply limit. A common choice of quadratic cost function would b e J = ∆ x T N Λ f ∆ x N | {z } terminal cost + N − 1 X n =0 ∆ x T n Λ∆ x n + ∆ υ T n Φ∆ υ n | {z } running cost (31) with ∆ x n := x n − x ? n state error b etw een actual states and reference state target x ? n , ∆ υ n := υ n − υ ? n con trol in- put error b etw een actual control signal and control signal target υ ? n corresp onding to reference state x ? n , horizon length N , and symmetric w eight matrices Λ f = Λ T f ≥ 0, Λ = Λ T ≥ 0, and Φ = Φ T > 0. Substituting (30) yields an equiv alent form of cost function as a quadratic func- tion of torques J (Υ) = c z + z Υ N + Υ T N Z Υ N (32) with c z := ( ¯ x 0 + ¯ g − X ? N ) T ˜ Λ( ¯ x 0 + ¯ g − X ? N ) + Υ ? T N ˜ ΦΥ ? N , z := 2( ¯ x 0 + ¯ g − X ? N ) T ˜ Λ ˆ H − 2Υ ? T N ˜ Φ , Z := ˆ H T ˜ Λ ˆ H + ˜ Φ , ˜ Λ := diag(Λ , · · · , Λ | {z } N − 1 copies , Λ f ) , ˜ Φ := diag(Φ , · · · , Φ | {z } N copies ) , Υ N :=      υ 0 . . . υ N − 1      , ¯ x 0 :=      F 1 . . . F N      x 0 , ¯ g := ˆ F      g . . . g      , X ? N :=      x ? 1 . . . x ? N      , Υ ? N :=      υ ? 0 . . . υ ? N − 1      , ˆ H := ˆ F        H H · · · H 0 H · · · H . . . . . . . . . . . . 0 0 · · · H        | {z } N × N copies , and ˆ F :=          I 0 · · · 0 F I . . . . . . . . . . . . 0 F N − 1 F N − 2 · · · I          . Similarly , one can see from the direct application of (30) that linear constraints on states in the form a neq x k ≤ b neq ∀ k = 1 , . . . , N with a neq ∈ R n c × 2( n a + n u ) , b neq ∈ R n c , and ≤ denoting the comp onent-wise relation translates to a linear con- strain t on torques as A neq Υ N ≤ B neq 9 with A neq := D neq ˆ H , B neq := ¯ b neq − D neq ( ¯ g + ¯ x 0 ) , ¯ b neq :=  b T neq · · · b T neq  T , and D neq := diag( a neq , · · · , a neq | {z } N − 1 copies ) . Therefore, the optimal con trol problem can b e formu- lated as min Υ N J (Υ) = c z + z Υ N + Υ T N Z Υ N s.t. A neq Υ N ≤ B neq ( ˙ q a ) T n υ n + υ T n diag( ¯ R ) υ n ≤ P max , 0 ≤ n ≤ N − 1 (33) where ( ˙ q a ) n := S T ˙ q n denotes the velocities of actuated join ts at time step n , and ¯ R :=  R 1 k 2 t , . . . , R n a k 2 t  T . Prop osition 5 Power supply limit defines a quadr atic c onstr aint on tor ques Υ . Pro of 2 Observe that the solution to velo cities takes the form ˙ q n = S T q F n x 0 + n − 1 X i =0 F n − 1 − i ( H υ i + g ) ! , wher e the op er ator S T q = [0 , I ] T sele cts the lower half of the states (i.e., velo cities). Plugging this into (33) gives a quadr atic ine quality Υ T n        0 . . . 0 C 0 . . . . . . . . . 0 . . . 0 C n − 2 C T 0 . . . C T n − 2 Ω        | {z } E Υ n +        0 . . . 0 β        T | {z } χ T Υ n ≤ P max (34) with C i = 1 2 H T ( F n − 1 − i ) T S q S , Ω = diag( ¯ R ) , and β = S T S T q F n x 0 + n − 1 X i =0 F n − 1 − i g ! . The ine quality (34) c an b e lifte d to sp ac e of ful l c ontr ol action se quenc e Υ N by p adding the matrix E and the ve ctor χ acting on Υ n with zer os in the expr ession (34). Henc e, it c an b e c onclude d that the p ower limit imp oses a quadr atic c onstr aint on Υ N . Corollary 6 The optimal c ontr ol pr oblem in (33) is a quadr atic al ly c onstr aine d quadr atic pr o gr am (QCQP) [41]. Theorem 7 QCQP in (33) is a non-c onvex pr oblem. Pro of 3 First note that the c ost function is c onvex sinc e Z > 0 . F urther, line ar c onstr aints pr eserve c onvexity. Ther efor e, we r estrict the pr o of to showing that p ower c onstr aint is not c onvex. The ine quality (34) c an b e iden- tifie d with the fol lowing quadr atic form : Υ T n E Υ n + χ T Υ n ≤ P max (35) This c onstr aint is not c onvex if and only if E is not p osi- tive semidefinite. R est of the pr o of is done sep ar ately for the fol lowing two c ases : i) ¯ R 6 = 0 ii) ¯ R = 0 , c orr esp ond- ing to ful l mo del of p ower limit including ele ctric al losses and r e duc e d mo del without ele ctric al losses, r esp e ctively. i) Supp ose that the the or em holds. This is e quivalent to E b eing p ositive semidefinite, fr om which Schur Complement defines the fol lowing c onditions : E ≥ 0 ⇔ Ω ≥ 0 and      C 0 . . . C n − 1      Ω      C 0 . . . C n − 1      T ≤ 0 . (36) It is obvious that the first ine quality Ω ≥ 0 on the right hand side automatic al ly r efutes the ine quality next to it which le ads to a c ontr adiction with the pr op osition E ≥ 0 . As a r esult, it is c onclude d that the p ower supply c onstr aint is non-c onvex when ¯ R 6 = 0 . ii) F or a dynamic al system (29), the matrix E c an b e tr e ate d as a c ontinuous function of ¯ R sinc e its entries dep end c ontinuously on R i s. Given this fact, it fol- lows fr om The or em 3.1.1 in [42] that al l eigenvalues of E ar e c ontinuous functions of ¯ R , i.e., λ i ( E ( ¯ R )) ∈ C 0 for al l i ∈ 1 , ..., n a wher e λ i ( E ) denotes the i th eigenvalue of E and C 0 sp ac e of c ontinuous func- tions. F urthermor e, it is known fr om the pr evious c ase that E is not p ositive semidefinite which implies min i λ i ( E ( ¯ R )) < 0 (37) for ¯ R 6 = 0 . Ther efor e, in or der to show that the p ower supply c onstr aint is non-c onvex stil l in this c ase, it is sufficient to verify (37) when ¯ R = 0 . Now, supp ose that the c ontr ary is true, i.e., min i λ i ( E (0)) ≥ 0 . (38) On the other hand, it fol lows fr om the fact tr(E( ¯ R)) = n × n a X i=1 λ i (E( ¯ R)) 10 that P i λ i ( E (0)) = 0 . Equation (38) and this r e- lation c an b e satisfie d simultane ously if and only if al l eigenvalues ar e zer o, which r e quir es E (0) to b e a nilp otent matrix, which is obviously not true, le ad- ing to a c ontr adiction with (38). Example 8 In or der to validate the algorithm pr op ose d in this se ction and demonstr ate its efficiency, a fin actu- ation system of a missile il lustr ate d in Fig. 5 is c onsid- er e d as an example system. The system c onsists of four identic al actuators e ach c onne cte d to a fin use d to c on- tr ol the missile by r e dir e cting aer o dynamic for c es on it. Ther efor e, an actuator of the fin i ∈ N a = { 1 , 2 , 3 , 4 } c ontr ols its p osition q i by applying tor ques u i which ar e subje ct to the fol lowing c onstr aints : (1) A ctuators have fixe d for c e limits (i.e., − ¯ u ≤ u i ≤ ¯ u asso ciate d with the p e ak curr ent r ating of their drive ele ctr onics. (2) A r e alistic tor que sp e e d curve is enfor c e d on actua- tors as a sp e e d dep endent for c e limit − u stall  1 + ˙ q i ˙ q max  ≤ u i ≤ u stall  1 − ˙ q i ˙ q max  wher e u stall and ˙ q max denote stal l for c e and no-lo ad sp e e d of actuators, r esp e ctively. (3) The system is supplie d by a c ommon p ower sour c e with maximum p ower output P max which defines a c onstr aint on tor ques dep ending on the p ower supply limit mo del use d by the c ontr ol lers. F or e ach fin actuation channel, we c onsider the me- chanic al system use d thr oughout Se c. 2 whose dynamics c an b e r epr esente d by the tr ansfer function (7). On the other hand, inter action of channels due to limits of the c ommon p ower supply c an b e mo dele d by formulating the entir e system as a multi DoF me chanic al system by cho osing q = [ q 1 q 2 q 3 q 4 ] T as gener alize d c o- or dinates, x = [ q T ˙ q T ] T as the system state, and υ = [ u 1 , u 2 , u 3 , u 4 ] T as a tor que input ve ctor, which yields a line ar system (28) in the state-sp ac e form with F c = " 0 4 × 4 I 4 × 4 0 4 × 4 ( − d/m ) · I 4 × 4 # , H c = " 0 4 × 4 (1 /m ) · I 4 × 4 # , and g c = 0 . Figure 5. Missile fin actuation system and its section view. T able 2 Numerical v alues of actuator and controller parameters. P arameter V alue m 1 kg m 2 d 0.05 Nm s / rad ¯ u 180 Nm u stall 500 Nm ˙ q max 4 rad / s ¯ R (0 . 0056) · [1 , 1 , 1 , 1] T Ohm A 2 / (Nm) 2 P max 750 W ∆ t 0.001 s N 300 Λ diag(2 , 2 , 2 , 2 , ∆ t, ∆ t, ∆ t, ∆ t ) · (0 . 5 / (∆ t ) 2 ) Φ diag(1 , 1 , 1 , 1) Λ f diag(0 . 1 , 0 . 1 , 0 . 1 , 0 . 1 , ∆ t, ∆ t, ∆ t, ∆ t ) · (10 / (∆ t ) 2 ) Thr e e optimal c ontr ol str ate gies ar e c onsider e d for this multi-DoF system. Contr ol lers adopt the same mo del for c onstr aints 1 and 2, i.e., an ine quality A neq Υ N ≤ b neq , wher e as they use differ ent mo dels for the c onstr aint 3 (i.e., p ower supply limit). In p articular, the fol lowing c ontr ol lers ar e c onsider e d : C1: The ar chite ctur e (33) unifying a c ontr ol ler and a dy- namic al lo c ation str ate gy is evaluate d first. T o that end, the limit on aggr e gate p ower c onsumption (i.e., Constr aint 2) is mo dele d by a quadr atic c onstr aint given in (34), henc e pr oviding a solution to b oth c on- tr ol and dynamic p ower al lo c ation pr oblems c oncur- r ently. C2: Mo difying the first c ontr ol ler by r e ducing the aggr e- gate p ower limit to p ower limit of individual actu- ators, a new c ontr ol ler c an b e obtaine d. This c on- tr ol ler c an b e formulate d as min Υ N J (Υ) = c z + z Υ N + Υ T N Z Υ N s.t. A neq Υ N ≤ B neq [ ˙ q n ] i [ υ n ] i + [ ¯ R ] i [ υ n ] 2 i ≤ ¯ P i , ∀ n ∧ ∀ i ∈ N a (39) wher e N a denotes the set of actuate d c o or di- nates, and [ ] i op er ator sele cts the i th actuator fr om this set. This c ontr ol ler r e quir es a pr e de- termine d static p ower al lo c ation solution to the e quation P 4 i =1 ¯ P i = P max . F or simulations, sinc e actuators ar e identic al, the trivial solution ¯ P i = P max / dim( N a ) ∀ i ∈ N a is use d. C3: The last c ontr ol ler c an b e derive d fr om the se c ond c ontr ol ler by mo deling an actuator’s p ower limit with the c onservative appr oximation (4). This enables formulation of the p ower limit in terms of a line ar c onstr aint, which c an b e c aptur e d by the optimal c on- tr ol ler min Υ N J (Υ) = c z + z Υ N + Υ T N Z Υ N s.t. A neq Υ N ≤ B neq [ υ n ] i ≤ ¯ P i / ˙ q max ,i , ∀ n ∧ ∀ i ∈ N a . (40) 11 -400 0 400 800 0 0.05 0.1 0.15 0.2 0.25 0 20 40 60 Figure 6. Cost-to-go function J (b ottom) scaled by a factor of 10 − 6 and total p ow er consumption P 4 i =1 P i of actuators (top), b oth as a function of time, for optimal controllers C 1 (solid-blue), C 2 (dashed-red), and C 3 (dotted-yello w). Using numeric al values given in T able 2 for the phys- ic al mo del and c ontr ol lers, simulations ar e p erforme d to evaluate c ontr ol lers for an example task of driv- ing the system state x fr om the initial c ondition x 0 = [0 . 5 , − 0 . 16 , 0 . 08 , 0 . 28 , 0 , 0 , 0 , 0] T to the tar get state x ? = 0 . Simulation r esults il lustr ate d in Figur es 6 and 7 indic ate fol lowings : • Contr ol ler C 3 , which adopts neither a dynamic al lo c a- tion str ate gy nor an exact p ower mo del, p erforms the worst, manifesting itself as the slowest settling time to tar get lo c ation and as the highest c ost function thr ough- out al l this time sinc e it c annot deliver mor e than half of the available p ower as shown in Fig. 6. This r esult is in favour of the main ide a pr op ose d in this p ap er, which is to exactly mo del p ower supply limit in c ontr ol lers. • Contr ol ler C 1 , i.e., the unifie d str ate gy pr op ose d in this se ction, outp erforms other str ate gies by pr oviding the minimum c ost-to-go function, as depicte d in Fig. 6. This is achieve d by al lo c ating lar ger p ortion of the available p ower to p articular joints which c ontribute to over al l tr acking p erformanc e mor e than other joints. This is in c ontr ast with the static al lo c ation which r e- serves e qual amount of p ower for e ach joint. In p ar- ticular, ther e ar e two me chanisms that c an act as a me dium for that str ate gy: 1) R esidual p ower of joints with lower priority ar e channel le d to joints with higher priority without violating the aggr e gate p ower limit so as to b enefit fr om the p ower supply as much as p ossible. This is il lustr ate d in Fig. 6 wher e c ontr ol ler C 1 uses 100% of the p ermissible p ower sour c e in the first 50 ms of the simulation. 2) L o c al p erformanc e of joints with lower priority is sacrific e d to impr ove over al l p erfor- manc e, me asur e d by the c ost-to-go function, as il lus- tr ate d in Fig. 7 wher e c ontr ol ler C 1 toler ates a slight delay in the r esp onse of actuator 3 to impr ove r esp onse of actuator 1. 4.2 Nonline ar Systems Ev en though dynamic allo cation increases control effi- ciency for linear systems as shown in the previous sec- tion, it requires solving a non-con vex QCQP optimiza- tion problem which is not computationally efficien t. F ur- thermore, if this finite-horizon optimization problem is generalized to nonlinear systems, this situation gets ev en w orse with non-quadratic and nonlinear ob jective and constrain t functions. In this subsection, an alternative approac h based on a control lyapuno v function (CLF) in conjunction with a p oint wise-minimum-norm con trol [43–45] is prop osed to obtain a more practical unified algorithm. Assuming that p ow er limit nonlinearity is satisfied, dy- namics in (9) can b e cast into state-space form ˙ x = f ( x ) + g ( x ) u (41) with x = " q ˙ q # , g ( x ) = " 0 M − 1 ( q ) # , and f ( x ) = " ˙ q − M − 1 ( q ) ( C ( q , ˙ q ) ˙ q + D ˙ q + G ( q )) # . No w, consider a task y ( x ) of joint p ositions with the goal to construct a feedback control law embedding a dynamic p ow er allo cation strategy and b eing compati- ble with constraints (e.g., torque limits and a p ow er sup- ply limit) suc h that dynamics of p ositional error ˜ y := y ( x ) − y ? ( t ) betw een the output function y ( x ) and a de- sired tra jectory y ? ( t ) is rendered GAS with an equilib- rium p oint e := h ˜ y , ˙ ˜ y i = 0. The prop osed control design is based on prior knowledge of a CLF which is defined as a contin uously-differentiable, prop er (i.e., radially- un b ounded), and p ositive-definite function V : R n → R + with n = dim( e ) such that inf u h ˙ V ( e ) i < 0 . (42) In other w ords, the existence of a CLF for trac king y ? ( t ) is equiv alent to the existence of a stabilizing controller whic h satisfies the inequality (42). Suc h a controller and a CLF pair can b e found b y feedback linearization [46] systematically as follo ws : Considering ˜ y ( x, t ) as the out- put function, input-output dynamics ha ve a relativ e de- gree of tw o, hence taking the form ¨ ˜ y = L 2 f ˜ y + ( L g L f ˜ y ) u − ¨ y ? (43) 12 -0.15 0 0.15 0.3 0.45 0 0.1 0.2 -180 -90 0 90 180 0 0.1 0.2 0 0.1 0.2 0 0.1 0.2 Figure 7. Position (top) and torque (bottom) tra jectories of actuators whose settling times are marked with circles for optimal con trollers C 1 (solid-blue), C 2 (dashed-red), and C 3 (dotted-yello w). Settling time is defined as the time required for the resp onse curv e to reach and stay within 5% of its initial v alue. where L X T denotes the Lie deriv ative of a tensor field T , whic h can b e a scalar function, a v ector field, or a one- form, along the flow of a vector field X [47]. Observe that error dynamics (43) is a time-v arying nonlinear system whic h can b e represented in state-space form d d t " ˜ y ˙ ˜ y # |{z} e = " ˙ ˜ y L 2 f ˜ y − ¨ y ? ( t ) # | {z } ˜ f ( t,e ) + " 0 n × 1 L g L f ˜ y # | {z } ˜ g ( t,e ) u. (44) As an illustrativ e example, consider the full joint-space con trol i.e., y ( x ) = q . In this case, input-output system (44) b ecomes a trivial extension of the original state- space system (41) with a time-dep endent term as ˜ f ( t, e ) = f ( x ) − h ( ˙ q ? ) T ( ¨ q ? ) T i T ˜ g ( t, e ) = g ( x ) , whic h sho ws that (44) is sufficien tly general to express a v ariet y of tasks. Therefore, a feedbac k-linearizing control la w u = ( L g L f y ( x )) − 1  ˜ u + ¨ y ? − L 2 f y ( x )  , (45) with an auxiliary con troller ˜ u transforms nonlinear dy- namics of the system (44) to those of a linear system d e d t = " 0 n × n I n × n 0 n × n 0 n × n # e + " 0 n × n I n × n # ˜ u. whic h can b e globally asymptotically stabilized by c ho osing auxiliary part of the controller as ˜ u = − K p ˜ y − K d ˙ ˜ y . (46) Resulting closed-lo op error dynamics can b e expressed in state-space form as d e d t = " 0 n × n I n × n − K p − K d # | {z } A cl e (47) A candidate CLF for the system in (44) is V ( e ) = e T P e where P is a symmetric positive-definite matrix solution to the Lyapuno v equation A T cl P + P A cl + W = 0 (48) for some symmetric p ositive-definite matrix W . Then, deriv ativ e of V along tra jectories of the closed-lo op sys- tem (47) satisfies ˙ V = − e T W e ≤ − εV ≤ 0 , (49) v erifying the existence of a exp onentially stabilizing CLF for the system (44) such that V ( e ( t )) ≤ V ( e (0)) e − εt with ε := λ min ( W ) /λ max ( P ). The abov e pro cedure pro vides a developmen t framework to build con trollers for nonlinear mechanical systems de- spite its w ell-known drawbac ks [43, 45]. In particular, a feedbac k linearization controller may apply unnecessar- ily large control signals by blindly cancelling nonlinear- ities that might b e b eneficial to conv ergence of the sys- tem to the desired tra jectory . Moreov er, that increase in con trol effort may conflict with ph ysical constraints 13 suc h as torque limits. Motiv ated from these disadv an- tages, [43, 44], and [48] prop osed a new family of con- trollers based on p oint wise solution to the QP problem min u ( u − u 0 ) T Φ( u − u 0 ) s.t. L ˜ f ( e T P e ) + L ˜ g ( e T P e ) u ≤ − e T W e (50) whic h seeks a minimum norm solution to control sig- nals around a baseline v alue u 0 w eigh ted b y a matrix Φ = Φ T > 0 under the stabilization constrain t expressed as an upp er bound on time-deriv ative of the CLF (49). The baseline control signal can b e c hosen in a v ariet y of wa ys. F or instance, u 0 = 0 can b e used to minimize the instantaneous control actions or u 0 = u ? to en- force a particular con trol p olicy where u ? is computed b y adding output of a baseline feedback law and feed- forw ard torques for trac king the desired tra jectory or ev en u 0 = u k − 1 to minimize the c hange in control ac- tions. Note that the QP problem is solved at each sam- pling time, and controller output is up dated accordingly . This strategy (hereinafter referred to as CLF-QP) has b een successfully implemented for locomotion con trol of bip edal robots [49] and for adaptive cruise con trol [50]. It can also serv e as a framework to formulate a unified algorithm combining controls and dynamic p ow er allo- cation by incorp orating physical constraints in to the op- timization framework (50). As explained in Subsection 4.1, there are tw o types of constraints : 1) Affine con- strain ts suc h as torque limits and speed/torque curve of eac h join t. These constraints can b e enforced by a linear inequalit y A neq ( q , ˙ q ) u ≤ b neq ( q , ˙ q ) that needs to b e sat- isfied at all times. 2) Instan taneous p ow er supply limit whic h can b e mo deled as u T Ω u + ˙ q T u ≤ P max . (51) With these constrain ts, the QP formulation of a unified algorithm takes the form min u ( u − u 0 ) T Φ( u − u 0 ) s.t. L ˜ f ( e T P e ) + L ˜ g ( e T P e ) u ≤ − e T W e A neq ( q , ˙ q ) u ≤ b neq ( q , ˙ q ) u T Ω u + ˙ q T u ≤ P max . (52) W e now state the main result of this subsection. Theorem 9 The unifie d algorithm b ase d on CLF-QP formulation (52) is a c onvex optimization pr oblem. Pro of 4 Assume that p osition and velo city me asur e- ments ar e available. The pr o of is done sep ar ately for the fol lowing two c ases : i) ¯ R 6 = 0 ii) ¯ R = 0 , c orr esp onding to ful l mo del of p ower limit including ele ctric al losses and r e duc e d mo del without ele ctric al losses, r esp e ctively. i) In this c ase, (52) is a QCQP pr oblem sinc e the p ower limit c onstr aint is quadr atic due to nonzer o Ω . The c onvexity fol lows fr om the fact that Φ is a symmetric p ositive-definite matrix as define d in (50) and that Ω is a diagonal matrix with p ositive entries qualifying it as p ositive-definite symmetric matrix. ii) In this c ase, (52) r e duc es to a standar d QP without a quadr atic c onstr aint sinc e the p ower supply limit b e c omes an affine c onstr aint with Ω = 0 due to ¯ R = 0 . Convexity is automatic al ly obtaine d as a r esult of this situation. Finally , observe that (52) might b e an infeasible QP problem since stabilization and ph ysical constraints may conflict with each other. T o ensure feasibilit y , we relax the QP problem by in tro ducing a slac k v ariable into CLF rate inequality such that intersection of feasible regions corresp onding to constraints is nev er empty . The relaxed v ersion of the unified algorithm can b e formulated as min u,p s ( u − u 0 ) T Φ( u − u 0 ) + c s p 2 s s.t. L ˜ f ( e T P e ) + L ˜ g ( e T P e ) u ≤ − e T W e + p s A neq ( q , ˙ q ) u ≤ b neq ( q , ˙ q ) u T Ω u + ˙ q T u ≤ P max . (53) with a slack v ariable p s and a p enalty co efficient c s as- so ciated with p s . Example 10 Consider the example in Se c. 3.2.2 wher e a two-link manipulator il lustr ate d in Fig. 3 is use d to evaluate a p assivity-b ase d c ontr ol ler. This time, the r ob ot example is use d to evaluate the pr op ose d CLF-QP-b ase d appr o ach that c ombines the dynamic al lo c ation and a task c ontr ol ler. T o that end, we assume that the r ob ot’s phys- ic al p ar ameters given in T able 1 ar e the same. As for the c ontr ol ler, in c ontr ary to that example, PD gains K p and K d of the fe e db ack-line arizing c ontr ol ler (45) ar e chosen as K p = ω 2 n K d = 2 ζ ω n , with ω n = 2 π (2 . 2) · I 2 × 2 and ζ = ( √ 3 / 2) · I 2 × 2 . This defines the close d-lo op system (47) whose CLF V = e T P e c an b e p ar ametrize d as P = " 2 ζ ω 2 n 2 ω n p 1 − ζ 2 2 ω n p 1 − ζ 2 2 ζ # . F urthermor e, ele ctric al losses of the motor and two phys- ic al c onstr aints ar e include d into the system mo del as fol- lows : (1) The motor r esistanc e normalize d by tor que c onstant is chosen as ¯ R = h 0 . 0833 , 0 . 222 i mOhm A 2 / (Nm) 2 . 14 (2) Joint tor que limits − ¯ u i ≤ u i ≤ ¯ u i ar e enfor c e d as        1 0 − 1 0 0 1 0 − 1        | {z } A neq " u 1 u 2 # | {z } u ≤        ¯ u 1 ¯ u 1 ¯ u 2 ¯ u 2        | {z } b neq with ¯ u = h 2000 , 1000 i Nm. (3) Power supply limit is enfor c e d as (51) with P max = 1 kW. Base d on this system mo del, thr e e c ontr ol lers ar e imple- mente d in a simulation envir onment for validation and p erformanc e c omp arison. C1: The first c ontr ol ler is CLF-QP-b ase d unifie d algo- rithm (53) pr op ose d in this subse ction with a zer o b aseline c ontr ol signal u 0 = h 0 , 0 i , identity weight- ing matrix Φ = I 2 × 2 and the p enalty c o efficient c s = 5 · 10 4 , r esp e ctively. C2: The se c ond c ontr ol ler is again a CLF-QP-b ase d uni- fie d algorithm but with static al lo c ation as opp ose d to the c ontr ol ler C 1 employing a dynamic al lo c a- tion. This differ enc e c an b e c aptur e d by mo difying QP pr oblem (53) as min u,p s ( u − u 0 ) T Φ( u − u 0 ) + c s p 2 s s.t. L ˜ f ( e T P e ) + L ˜ g ( e T P e ) u ≤ − e T W e + p s A neq ( q , ˙ q ) u ≤ b neq ( q , ˙ q ) ¯ R 1 u 2 1 + ˙ q 1 u 1 ≤ P max / 2 ¯ R 2 u 2 2 + ˙ q 2 u 2 ≤ P max / 2 . The CLF and other p ar ameters of this c ontr ol ler ar e the same as those use d by c ontr ol ler C 1 . C3: The thir d str ate gy is the standar d fe e db ack- line arization c ontr ol ler describ e d in this se ction. Using initial c onditions x (0) = h − π / 2 0 0 0 i , simula- tions ar e c onducte d to evaluate these thr e e c ontr ol lers in a joint-sp ac e r e gulation task towar d a tar get p osition q ? = h π / 2 0 i . This c orr esp onds to the output function y ( x ) = q − q ? which further defines the c ontr ol pr o c e dur e explaine d in this se ction. Simulation r esults il lustr ate d in Figur es 9 and 8 show the fol lowings : Contr ol ler C 1 demonstr ates the b est p erformanc e. This c an b e se en fr om the p osition tr aje ctories of joints as wel l as the CLF. In p articular, joint 1 r e aches the tar get angle in the shortest time while pr oviding the minimum amount of deviation for joint 2 fr om its desir e d p osition. This is achieve d by channel ling gr e ater amount of p ower to joint 1 than other c ontr ol lers without exc e e ding tor que limits. On the other hand, static al lo c ation str ate gies ar e c omp ar able to e ach other. This manifests itself as CLF tr ends with similar time c onstants. However, the c ontr ol ler C 2 pr o duc es a less deviation in joint 2 wher e as the c ontr ol ler C 3 pr o- vides a shorter settling time for joint 1. While doing so, the c ontr ol ler C 2 r e quir es less c ontr ol effort b e c ause of the fact that it pr ovides the same ne gative CLF r ate with the the minimum tor que. 0 1 2 3 0 0.5 1 -4 -2 0 2 Figure 8. V alue of CLF V (top) scaled by a factor of 10 − 3 and total p ow er consumption P 2 i =1 P i of joints (b ottom), b oth as a function of time, for controllers C 1 (solid-blue), C 2 (dashed-red), and C 3 (dotted-purple). -1 0 1 -2 0 2 0 0.5 1 -6 -3 0 3 0 0.5 1 Figure 9. P osition (top), torque (middle), and p ow er (bot- tom) tra jectories of joint 1(left) and joint 2 (right) for con trollers C 1 (solid-blue), C 2 (dashed-red), and C 3 (dot- ted-purple). 15 5 Preliminary Experimental Results In this section, w e present results of preliminary ex- p erimen ts conducted on a single DoF electromec han- ical actuator that we built to comparativ ely ev aluate t w o controllers with exact and approximate models of p o w er supply limit, resp ectively , as detailed in Sec. 3.2. The exp erimental setup illustrated in Fig. 10 includes the actuator, a p ow er-supply , and a test rig, which is a self-contained system capable of applying controlled external loads to and measuring output p osition and p o w er consumption of the actuator. The actuator con- sists of a drivetrain, a PMSM with an encoder, and con- trol/driv er electronics based on an FPGA with a generic six bridge p o w er amplifier. The electronics run a field orien ted controller for torque control and a PID-type tra jectory tracking controller with CI that is designed according to the procedure in Sec. 3.2. Parameters of the actuator and other details of the exp erimental setup are giv en in T able 3. Tw o con trollers are implemented on the actuator for exp erimen tal ev aluation. While b oth controllers are of the PID-t yp e, they use different mo dels of the p ow er supply limit nonlinearity , as mentioned abov e. The first algorithm, which is considered as a baseline controller and abbreviated as C 1 from now on, employs the exact mo del (3), which ignores the electrical losses of the mo- tor since its efficiency is close to unit y , in combination with the torque saturation nonlinearity corresp onding to the peak curren t of the driver to find a time-v arying instan taneous torque limit as u max ( t ) = min( I max k t , P max / ˙ q ( t )) . On the other hand, the second controller which will b e referred as C 2 from no w on, uses the appro ximate mo del (4), which defines a standard torque saturation sat( u, u max ) with u max = min( I max k t , P max / ˙ q max ) = P max / ˙ q max = 100 N m. T able 3 P arameters of the actuator, controllers, and the p ow er supply in the exp erimental setup. P arameter V alue Motor Maxon EC-4p ole 30 Phase-to-Phase Resistance 0.1 Ω Motor Efficiency 95 % V oltage 24 V Inertia 1 kg m 2 Viscous damping 0.05 N m / (rad/s) T orque control sampling rate 20 kHz T orque control bandwidth 500 Hz Maxim um sp eed ˙ q max 4 rad/s P eak driv er curren t I max 32 A T orque constant k t 6 N m / A P o w er supply limit P max 400 W P osition con trol sampling rate 2 kHz Figure 10. Exp erimental setup. Tw o t ypes of exp erimen ts are conducted to v alidate con- trollers and to compare their resulting closed-lo op con- trol p erformance : 1) Time-domain exp eriments : Step-input in p osition- command with differen t amplitudes are applied to con- trollers whose time-domain resp onses are illustrated in Fig. 11 from which a comparison can b e made in terms of their maximum o v ersho ots and settling times. In par- ticular, as is evident from the results rep orted in T a- ble 4, the baseline controller C 1, which uses the exact mo del of p ow er limitation, outp erforms the controller with the appro ximate model clearly b y pro viding nearly half settling-time for the command with the highest am- plitude without an y significan t increase in the o v ersho ot whereas resp onses of controllers b ecome less differen t as amplitude decreases. This is due to the fact that larger torques are needed to achiev e the desired resp onse to step commands with greater amplitudes, which con- troller C 2 cannot deliver because of the conserv ative lim- itations on admissible torques defined b y the appro xi- mate p ow er limit mo del. 0 2 4 -200 0 200 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 -100 0 100 200 Figure 11. Po wer consumptions (bottom), motor torques (middle) and p osition responses (top) of controllers C 1 (solid blue line) and C 2 (dashed red line) to 1 ◦ ,2 ◦ , and 3 ◦ step commands (dash dotted green line in the top figure). T orque limits due to the driv er’s p eak curren t rating and pow er sup- ply limit are marked with the dotted orange lines. 16 T able 4 Settling time and P ercen t Overshoot (PO) v alues of con- trollers for different step commands. Con troller Amplitude Settling Time PO C 1 1 ◦ 0.036 0 2 ◦ 0.039 0 3 ◦ 0.043 0 C 2 1 ◦ 0.038 0 2 ◦ 0.051 6.0 3 ◦ 0.071 24.30 2) F requency-domain exp eriments: W e compare the closed-lo op frequency resp onse of controllers by stim - ulating the system with a chirp command which has a unit amplitude and frequency conten t in the range b et w een 1 H z and 24 H z and by measuring the result- ing frequency resp onse functions, which are illustrated in Figure 12. The Bo de plots show that controller C 1 yields a b etter closed-lo op tracking performance by pro viding higher magnitude and less phase dela y in the high frequency region (esp ecially ab o v e 10 Hz). F ur- thermore, the p ow er consumption of b oth controllers giv en in Figure 13 shows that the controller C 1 chan- nels as m uch pow er as p ermissible to the motor without exceeding the p ow er supply limit by reacting quickly to the changes in the measurements. This result v ali- dates not only the effectiv eness of the prop osed control approac h but also its safety . -10 -5 0 Figure 12. Magnitude (top) and phase (b ottom) resp onse of the actuator with con trollers C 1 (solid blue) and C 2 (dashed orange). 0 10 20 -200 0 200 0 10 20 Figure 13. Po wer consumption of con trollers C 1 (left) and C 2 (right) under pow er supply limit (dotted yello w) during the chirp test. In conclusion, b oth time-domain and frequency-domain exp erimen ts on a simple 1-DoF actuator provide pre- liminary evidence verifying theoretical adv antages of the main idea in this pap er, whic h is to use exact mo del of p o w er supply limit in con trollers instead of the appro x- imate mo del. 6 Conclusions and F uture W orks In this work, motion con trol systems with strict pow er limits are considered. Effects of these physical limits as b eing a nonlinearity on dynamics are analysed rigorously in b oth frequency and time domains. In the former anal- ysis, our exploration of the frequency response b ehavior of p ow er supply limits follows from deriv ations of de- scribing function and maximum bandwidth frequency . Time domain inv estigation is done with reference to con- troller design. In particular, we consider dynamic clas- sical con trollers for linear and nonlinear Euler-Lagrange systems. Finally , t w o general c on trol arc hitectures are prop osed to improv e control p erformance b y leveraging the p ow er supply limit through dynamic p ow er alloca- tion for m ulti DoF systems. In particular, these archi- tectures are developed in the framew ork of constrained optimal con trol and constrained in v erse optimal con trol where QP problems are formulated for linear systems and nonlinear systems, resp ectively . Finally , preliminary results of an experiment are presented and briefly dis- cussed to supp ort the theoretical contributions in the pap er. Since the presented material provides a new p ersp ec- tiv e on an old problem which has been b ecoming more imp ortan t with recent developmen ts in mobile systems, there are v arious future work directions : Analysis on stabilit y and p erformance limitations due to p ow er limit can b e geared tow ards underactuated nonlinear systems for which swing-up energy control is a com- monly used strategy [51, 52]. The connection b etw een passivit y-based control and p o w er supply limits might b e in teresting to explore further since passivity-based con trol is related to regulation of system energy . Con- v exification of the unified controller which solves b oth dynamic p ow er allo cation and con trol problems for linear mechanical systems b y a finite-horizon QP prob- lem in conjuction with p ow er supply limit constraints is also an interesting av enue to explore since it would expand the application areas of those controllers by reducing computational requirements. Finally , further exp erimen tal results are needed to verify b oth classi- cal controllers considering pow er supply limits and the unified algorithm em b edding dynamic p o wer allo cation and feedback control. References [1] I. Kelly , O. Holland, C. Melh uish, SlugBot: A Robotic Predator in the Natural W orld, in: Proc. of the 5th Int. 17 Symposium on Artificial Life and Rob otics, 2000. [2] G. A. Landis, Exploring Mars with Solar-Po wered Rov ers, in: Conf. Record of the 31st IEEE Photo voltaic Sp ecialists Conf., 2005. [3] R. W elch, D. Limonadi, R. Manning, Systems Engineering the Curiosity Rov er: A Retrosp ective, in: Proc. of the 8th Int. Conf. on System of Systems Engineering, 2013. [4] H. Kolven bach, M. Hutter, Life Extension: An Autonomous Docking Station for Recharging Quadrup edal Rob ots, in: Proc. of the 11th Int. Conf. on Field and Service Robotics, 2017. [5] G. D. Kenneally , A. De, D. E. Ko ditschek, Design Principles for a F amily of Direct-Drive Legged Rob ots, IEEE Rob otics and Automation Letters 1 (2) (2016) 900–907. [6] A. Abate, J. W. Hurst, R. L. Hatton, Mechanical Antagonism in Legged Rob ots, in: Pro c. of the Rob otics : Science and Systems, 2016. [7] S. Rezazadeh, A. Abate, R. L. Hatton, J. W. Hurst, Rob ot Leg Design : A Constructive F ramework, IEEE Access 6 (2018) 54369 – 54387. [8] S. Collins, A. Ruina, R. T edrake, M. Wisse, Efficient Bipedal Rob ots Based on Passiv e-Dynamic W alkers, Science 307 (5712) (2005) 1082–1085. [9] G. Secer, U. Saranli, Control of Planar Spring-Mass Running Through Virtual T uning of Radial Leg Damping, IEEE T rans. on Rob otics 34 (5) (2018) 1370–1383. [10] S. Seok, A. W ang, D. Otten, S. Kim, Actuator Design for High F orce Proprio ceptive Control in F ast Legged Lo comotion, in: Proc. of the IEEE/RSJ Int. Conf. on Intelligen t Rob ots and Systems, 2012. [11] J. C. Doyle, R. S. Smith, D. F. Enns, Control of Plants with Input Saturation Nonlinearities, in: Pro c. of the American Control Conf., 1987. [12] J. W. Choi, S. C. Lee, Antiwindup Strategy for PI-Type Speed Con troller, IEEE T rans. on Industrial Electronics 56 (6) (2009) 2039–2046. [13] H. A. F ertik, C. W. Ross, Direct Digital Control Algoritm With Anti-Windup F eature, ISA T ransactions 6 (4) (1967) 317–328. [14] R. Hanus, M. Kinnaert, J. L. Henrotte, Conditioning T echnique : A General Anti-Windup and Bumpless T ransfer Method, Automatica 23 (1987) 729–739. [15] E. Gilb ert, I. Kolmanovsky , Nonlinear T racking Con trol in the Presence of State and Control Constraints : A Generalized Reference Governor, Automatica 38 (12) (2002) 2063–2073. [16] I. Kolmanovsky , E. Garone, S. D. Cairano, Reference and Command Governors : A T utorial on Their Theory and Automative Applications, in: Pro c. of the American Control Conf., 2014. [17] A. S. Hodel, C. E. Hall, V ariable-Structure PID Control to Preven t Integrator Windup, IEEE T rans. on Industrial Electronics 48 (2) (2001) 442–451. [18] S. Galeani, S. T arb ouriech, M. T urner, L. Zaccarian, A T utorial on Mo dern Anti-Windup Design, Europ ean Journal of Control 3 (4) (2009) 418–440. [19] K. J. Astrom, L. Rundqwist, Integrator Windup and How to Avoid It?, in: Proc. of the American Control Conf., 1989. [20] G. Grimm, J. Hatfield, I. Postleth waite, A. R. T eel, M. T urner, L. Zaccarian, Anti-Windup for Stable Linear Systems with Input Saturation : An LMI-Based Synthesis, IEEE T rans. on Automatic Control 48 (9) (2003) 1509–1525. [21] E. F. Mulder, M. V . Kothare, M. Morari, Multiv ariable Anti-Windup Controller Synthesis Using Linear Matrix Inequalities, Automatica 37 (9) (2001) 1407–1416. [22] C. Pittet, S. T arb ouriech, C. Burgat, Stability Regions for Linear Systems with Saturating Controls via Circle and Popov Criteria, in: Proc. of the 36th Conf. on Decision and Control, 1997. [23] H. Hindi, S. Boyd, Analysis of Linear Systems with Saturation using Convex Optimization, in: Pro c. of the 37th IEEE Conf. on Decision and Control, 1998. [24] D. W u, R. T o dd, A. J. F orsyth, Adaptive Rate-Limit Con trol for Energy Storage Systems, IEEE T rans. on Industrial Electronics 62 (7) (2015) 4231–4240. [25] A. Gelb, W. E. V. V elde, Multiple-Input Describing F unctions and Nonlinear System Design, McGraw-Hill, 1968. [26] H. Li, J. Shang, B. Zhang, X. Zhao, N. T an, C. Liu, Stability Analysis With Considering the T ransition Interv al for PWM DC-DC Conv erters Based on Describing F unction Metho d, IEEE Access 6 (2018) 48113–48124. [27] I. Boik o, Analysis of Sliding Mo des in the F requency Domain, International Journal of Control 50 (9) (2009) 1442–1446. [28] S. F. Hsiao, D. Chen, C. J. Chen, H. S. Nien, A New Multiplefrequency Small-Signal Mo del for High-Bandwidth Computer V-core Regulator Applications, IEEE T ransactions on Pow er Electronics 31 (1) (2016) 733–742. [29] U. P erez-V entura, L. F ridman, Design of Sup er-Twisting Control Gains : A Describing F unction Based Metho dology , Automatica 99 (2019) 175–180. [30] J. Munkres, Section 26 in T op ology, 2nd Edition, P earson, 2000. [31] M. Reyhanoglu, A. v an der Schaft, N. Mcclamroch, I. Kolmanovsky , Dynamics and Control of a Class of Underactuated Mechanical Systems, IEEE T rans. on Automatic Con trol 44 (9) (1999) 1663–1671. doi:10.1109/ 9.788533 . [32] R. Ortega, A. Loria, P . J. Nicklasson, H. Sira-Ramirez, Passivit y-Based Control of Euler-Lagrange Systems : Mechanical, Electrical, and Electromec hanical Applications, Springer-V erlag, 1998. [33] R. M. Murray , Z. Li, S. S. Sastry , Section 5.4 in A Mathematical In troduction to Rob otic Manipulation, CRC Press, 1994. [34] A. Visioli, Section 3.3.2 in Practical PID Control, Springer, 2006. [35] L. Zaccarian, A. R. T eel, Section 2.3.5 in Mo dern Anti- windup Synthesis, Princeton U niversit y Press, 2011. [36] T. Hu, Z. Lin, B. M. Chen, An Analysis and Design Method for Linear Systems Sub ject to Actuator Saturation and Disturbance, Automatica 38 (2) (2002) 351–359. doi: 10.1016/s0005- 1098(01)00209- 6 . [37] C. Lin, Q. W ang, T. H. Lee, An Improv ement on Multiv ariable PID Controller Design via Iterative LMI Approach, Automatica 40 (3) (2004) 519–525. [38] D. E. Ko ditschek, The Application of T otal Energy as a Lyapuno v F unction for Mechanical Control Systems, T ech. rep., University of Pennsylv ania (1989). [39] P . T omei, Adaptive PD Con troller for Robot Manipulators, IEEE T rans.on Rob otics and Automation 7 (4) (1991) 565– 570. [40] H. K. Khalil, Section 4.2. Nonlinear Systems, Prentice Hall, 2002. 18 [41] S. Boyd, L. V andenberghe, Section 4.4 in Con vex Optimization, Cambridge Universit y Press, 2004. [42] J. M. Ortega, Numerical Analysis : A Second Approach, V ol. 3 of Classics in Applied Mathematics, SIAM, 1990. [43] R. F reeman, P . Kokoto vic, Optimal Nonlinear Con trollers for F eedback Linearizable Systems, in: Proc. of the American Control Conference, 1995. [44] R. A. F reeman, P . V. Kokoto vic, Inv erse Optimalit y in Robust Stabilization, SIAM Journal on Control and Optimization 34 (4) (1996) 1365–1391. [45] A. D. Ames, K. Gallowa y , J. W. Grizzle, K. Sreenath, Rapidly exp onentially stabilizing control lyapuno v functions and hybrid zero dynamics, IEEE T ransactions on Automatic Control 59 (4) (2014) 876–891. [46] K. Y. Chao, M. J. Po well, A. D. Ames, P . Hur, Unification of Lo comotion Pattern Generation and Control Lyapuno v F unction-Based Quadratic Programs, in: Proc. of the American Control Conference, 2016. [47] J. M. Lee, Chapter 9 in Introduction to Smooth Manifolds, Springer, 2012. [48] J. A. Primbs, V. Nevistic, J. C. Doyle, Nonlinear Optimal Control : A Control Lyapunov F unction and Receding Horizon P erspective, Asian Journal of Con trol 1 (1) (1999) 14–24. [49] K. Gallow ay , K. Sreenath, A. D. Ames, J. W. Grizzle, T orque Saturation in Bip edal Rob otic W alking Through Control Lyapuno v F unction-Based Quadratic Programs, IEEE Access 3 (2015) 323–332. [50] A. Mehra, W. Ma, F. Berg, P . T abuada, J. W. Griz zle, A. D. Ames, Adaptive Cruise Control: Exp erimental V alidation of Adv anced Con trollers on Scale-Mo del Cars, in: Pro c. of the American Control Conf., 2015. [51] M. W. Spong, The Swing Up Control for the Acrob ot, IEEE Control Systems Magazine 15 (1) (1995) 49–55. [52] K. J. Astrom, K. F uruta, Swinging up a p endulum by energy control, in: Proc. of the IF AC 13th W orld Congress, 1996. 19