Nonlinear Trajectory Optimization Models for Energy-Sharing UAV-UGV Systems with Multiple Task Locations

Nonlinear T rajectory Optimization Models f or Energy-Sharing U A V -UGV Systems with Multiple T ask Locations Minsen Y u a n, Amanuel Adane, James Humann, and Y ue Y u Abstract — Energy-sharing U A V -UGV systems extend the en- durance of Uncrewed Aer ial V ehicles (U A Vs) by levera ging Uncrewed Ground V eh icles (UGVs) as mo bile char ging sta- tions, enablin g persistent autonomy in in f rastructure-sparse en vironments. T rajectory optimization for these systems is often challenging due to UGVs’ terrain access constraints and the dis- crete nature of task scheduling. W e p ropose a smooth nonlin ear program model f or th e joint trajectory optimization f or these systems. Unlike existing models, the pr oposed m odel allo ws smooth parameterization of UGVs’ terrain access constrain t s and supports partial U A V rechar ging. Further , it introduces a smooth app roxima tion of di sjunctive constraints that el i minates the need for comput ationally expensive integer programming and enables efficien t solutions via nonlin ear pr ogramming algorithms. W e demo nstrate the proposed m odel on a one- U A V -one-UGV system with multiple task locations. Compared with mixed-integer nonlin ear programs, this model re duces the computation time by orders of magnitude. I . I N T R O D U C T I O N Collaborative UA V -UGV systems p rovide a powerful plat- form to integrate the heterogeneo u s strength s of aerial and groun d robots v ia ene rgy-sharing . In these systems, Un- crewed Ground V ehicles (UGVs) serve as m obile do cking and charging stations that ef fectiv ely extend e d the endur ance of Uncr ewed Aerial V eh icles (U A Vs) via sharing energy stored on UGVs. By comb ining UGVs’ superior paylo a d ca- pacity together with U A Vs’ agility and elev ated per spectiv e, these system s en able a wide r ange of robotic applications in remo te and challengin g en viron ments where access to traditional energy infrastructure is sparse. Examples of these applications inclu d e mapping, disaster resp onse, surveillance, and precision agricultur e [1]–[3]. Coordinatin g the trajectories for energy-shar in g U A V - UGV systems presents a un ique challenge, especially in missions with multiple tasks loca tions. This difficulty is twofold. First, The UGV’ s terrain access constraints restrict the UGVs to a k n own road network. Second, sched u ling the the order in which the U A Vs and UGVs completes th eir individual tasks often requ ire optimizin g discrete decisions. When these discrete decisions ar e co upled with the ter rain access con straints, the resulting optimization is often non- conv ex an d co mputation a lly expensive to solve. One approach to o p timize UA V -UGV trajectories is via the two-echelon r outing p r oblem , wh e re the traje c to ries of M. Y ua n and Y . Y u are with the Depa rtment of Ae rospace Engineering and Mechanics, Uni ve rsity of Minnesota, Minneapolis, MN 55455, USA ( { yuan0450, yuey } @umn.edu ). A. Adane is with the Duffield Colle ge of Engineering, Corne ll Uni versit y , Itha ca, NY 14853, USA ( aa2888@cornell .edu ). J. Humann is with DEVCOM Army Re- search Laboratory ( james.d.huma nn.civ@army.mil ). Y . Y u would like to thank Samet Uzu n for he lpful early discussions. the U A V and UGV are computed in two distinct p hases. One appro a ch is to first determin ing the U A V trajectory , then o ptimize th e UGV trajectory to support the U A V [4 ]. Another approach is to first determine the UGV trajectory that provide a set of feasible r endezvous locations to support U A V r echarging, then p lan the U A V trajecto ry based these locations [5], [ 6]. Early results uses a gr e edy alg orithm to generate the rende z vous loc ations, then compute the U A V trajectory accor dingly [5], [7]. More recent results augment this approach with an async hronou s team framew ork [8] an d task allocation heuristics b ased o n m inimum set covering [6 ]. Another appr oach is to d ecouple the plan ning of UA V and UGV trajectory using reachable sets constraints [9]. In all such cases, the resulting trajectory quality depends o n how to partition the the plann ing between the U A V and UGV . T o dir ectly address the coupling b etween UA Vs and UGVs, an alternative app roach is to jointly o ptimize their trajectories. On approach for such optimization is mixed integer linear prog rams (MILP), which comp utes an optimal U A V -UGV trajectory via op timization over in teger variables [10], [11]. Since MILP is often comp utationally expe nsi ve, there have been se veral strategies to simplify the problem under stro n ger assum p tions. For examp le, if the potential rendezvous lo cations are known, genetic algorithm provides an efficient a lter nativ e to MILP [ 1 2]. If the order of U A V and UGV actions are kn own, then optimizing U A V -UG V trajectories red uces to a con vex seco nd-ord er cone progr am [13]. A mo re recent direction is deep reinf o rcement learning methods to train en coder-decoder based tran sformer net- works that generates the U A V -UGV trajector ies [14]–[ 16]. These trajectories can no t o nly minimize mission time but also account for stochastic battery usage by co nstraining on the risk of mid-m ission batter y depletion [17]. There are se veral limitatio ns in the existing mode ls for U A V–UGV trajec tory optimiz a tion. First, the UGV road network is often modeled as a g raph with discretized n odes, causing the plann ing com plexity to grow rapidly with the number of nod e s [4]–[ 8], [10]–[ 12]. Second, most existing models assume that the U A V battery is fully recha rged befo re takeoff [6 ] –[8], [12], [13], [15]. While this assum ption is compatible with battery replacement, it is not suitab le fo r wireless charging that allow partial recharge. Finally , MILP methods rely on integer variables to model discrete decisions, leading to solu tion times that g r ow expon entially with the number of integer v ariables [4] – [8], [10], [11]. T o address these limitations, we p ropose a novel n onlinear trajectory o ptimization mo del for energy-sharing UA V–UGV systems with m ultiple task location s. W e base this mo del on three key featu res. First, we m odel the UGV road network as the un ion of multiple continuo us p a ths joined at one junction point. This app roach avoids the need for nod e discretization along each p ath and allows smooth parame ter ization o f UGV’ s terrain access constraints. Second , we model the d y- namics of U A V battery via distinct charging and discharging curves. Compared with traditional mode ls that alw ays on ly allow fu ll recharges, this ap proach a llows p artial recharging and increase the flexibility of U A V trajec tories. T h ird, we de- velop a smo othing approach to model disjuncti ve co nstraints that are necessary for scheduling optimal task c o mpletion and battery rechargin g. W e first transform the disjunctiv e constraints into nonsmooth constraints defined by pointwise minimum f u nction, then propose a smooth ap proximatio n of these co nstraints based on ℓ p -norm s. This ap proxima tion enables efficient solutio n with smooth nonline a r optimizatio n methods, byp assing the need for comp utationally e xpen - si ve integer progr amming. W e demonstrate this mod el on a one-UA V–one-UGV system with multiple task locations. Compared with mixed-integer nonlin ear pr ogramm in g, this model reduces so lution time from hours to minu tes while maintaining robust success rates in numerical experiments. I I . T R A J E C T O RY O P T I M I Z AT I O N M O D E L S F O R E N E R G Y - S H A R I N G UA V - U G V S Y S T E M S W e m o del the trajec tory o ptimization fo r energy-sharing U A V -UGV system—where the UGV serve as mo bile ch arg- ing stations for the U A V —as a constrained optimization problem . This model includes the optimization variables used to p a r ameterize the U A V -UGV trajector y , the objective function , and the physical and operational co nstraints of th e U A V -UGV system. A. T r ajectory V ariables W e con sider three classes o f variables when o ptimizing th e U A V -UGV trajectory: variables shared by b oth trajectories, variables specific to the U A V trajector y , and variables specific to th e UGV trajec to ry . T o simplify th e notation, we let J k K : = { 1 , 2 , . . . , k } for any k ∈ N . a) Shar ed V ariables: Let N ∈ N den o te the total number of time stamps along the UA V trajectory , wh ich is the same as th e total nu mber o f time steps along the UGV trajectory . W e let s k ∈ R ≥ 0 denote the time duration between the k -th and the ( k + 1 ) -th time stamp for all k ∈ J N K . b) U A V V ariables: W e parameterize the U A V trajecto ry using its p lanar po sition and its remaining time of flight supported by the curr ent b attery le vel. W e let r A k ∈ R 2 denote the projection of UA V’ s 3D position onto th e xy-plane at th e k -th time stamp for all k ∈ J N K , and e k ∈ R ≥ 0 denote the remaining time of flight of the UA V at time step k . c) UGV V ariables: W e assume that th e UGV only moves on a road network in the shape of a star graph, consisting of m G arms that intersect on ly at a single ju n ction point. W e parame te r ize the UGV position trajectory using two variables. W e let r G k ∈ R 2 denote the planar position of the UGV at the k - th tim e stamp. In addition , we let p k ∈ R m G ≥ 0 denote the UGV’ s position o n the star graph with m G arms. Each entry of p k correspo n ds to one arm of the star graph, and the value of a positi ve entry is the distance the UGV trav eled f rom the graph center along that arm to its cu r rent position. At any time stamp, p k has at most one po siti ve entry , in dicating that the UGV occupies at most one a rm of the grap h at a time. B. Objective F unction W e choose the total mission time as the objective func- tion. Minimizing total tim e directly captures the oper ational efficiency of the UA V–UGV system, encourag in g timely completion of all required tasks wh ile im p licitly balan cing trav el, waiting, and charging time . Since th e time durations between co n secutiv e time stamps are optimization variables, the total mission time is g iven by P N − 1 k =1 s k . (1) C. T rajectory Co n straints W e co nsider three classes of constrain ts for the joint U A V–UGV trajectory: co nstraints th a t d epend o n ly o n U A V variables, con stra in ts that depend only o n UGV variables, and constraints that coup le bo th U A V and UGV variables. a) U A V Constr aints: W e let r 0 ∈ R 2 and r f ∈ R 2 denote the initial and final position of the U A V trajectory , respectively . W e con sid er the f ollowing initial an d final constraints r A 1 = r 0 , r A N = r f . (2) In add ition, we consid e r th e following con straint on the speed of the U A V   r A k +1 − r A k   2 ≤ v A max s k , (3) for all k ∈ J N − 1 K , where v A max ∈ R > 0 is the ma x - imum U A V speed. Note that the U A V po sition is three - dimensiona l, whereas the above constraint accoun ts only for planar motion. Th is is because we assume that the UA V operates at a con stan t altitude, and that the time req uired for vertical motion ( e.g. , durin g takeoff a n d landin g) is negligib le compare d with the time spent in ho rizontal motion. W e consider the case wh e re th e UA V must visit a set of task locations along its trajecto ry ( e .g., to monitor areas of interest). Let m A ∈ N deno te the total numb er of task locations, and let a 1 , . . . , a m A ∈ R 2 denote the locations of these task locatio ns on the xy -plan e. W e consider the following logical constraints: ∃ k ∈ J N K s.t. r A k = a i , (4) for all i ∈ J m A K . These constrain ts ensure that the U A V visits each task point at least on ce along its trajectory . b) UGV Con straints: W e consider the case where th e initial and final position of the UGV coincide with those of the U A V , namely , r G 1 = r 0 , r G N = r f . (5) The UGV can only move along a star-shaped road ne twork . At the k -th time stamp, r G k ∈ R 2 denotes the planar p osition of the UGV , and p k ∈ R m G ≥ 0 denotes its position along the arms of the star graph. I n particular, p k = 0 m G indicates that the UGV is lo c ated at the jun ction poin t (i.e., the center of the star g raph). Moreover , the UGV lies on the j -th a rm if and on ly if [ p k ] j > 0 . W e introd uce a non linear mapping g : R m G ≥ 0 → R 2 such that r G k = g ( p k ) . (6) Let p max ∈ R m G ≥ 0 denote a parameter vector whose j -th entry represents th e maximu m allowable distance of th e UGV from the junction point along the j - th arm, correspond ing to the furthest task p oint o n that arm. W e consider the following constraints on p k : 0 m G ≤ p k ≤ p max , r G k = g ( p k ) , p ⊤ k  1 m G 1 ⊤ m G − I m G  p k = 0 , (7) for all k ∈ J N K . These constraints ensure that the UGV occupies at most one arm of the star g r aph at any time stamp, and that its p lanar p osition r G k is coupled with th e vector p k throug h the nonlinear mapping g . Furthermo re, we c onsider the following constraints on the speed of the UGV k p k +1 − p k k 1 ≤ v G max s k (8) for all k ∈ J N − 1 K , where v G max ∈ R > 0 is the maximum UGV speed. The se con straint in (8), when combin e d with the co nstraint in (7), ensure tha t the distance tra veled by the UGV is upper bound by the p roduct of its maximum speed and the time used. Similar to the UA V , the UGV also must visit a set o f task locations. Since the UGV can only move on a roa d network modeled a s a star grap h, it is suf ficient to specify the furthest point along e a ch arm that the UGV m ust r each. Without lo ss of gen erality , we take the endpoint of each ar m to b e the correspo n ding furthest point. T o th is end, we consider the following constraints: ∃ k ∈ J N K s.t. b ⊤ j ( p k − p max ) = 0 , (9) for all j ∈ J m G K , wh ere b j ∈ R m G denotes the j -th colu m n of the id entity m a tr ix I m G . These constraints ensure that the UGV trajector y reaches the en d of each arm of the road network at least once along its trajectory . c) C oup ling Constraints: W e consider the following bound s on the time duration between two c o nsecutive time stamps, which constrain both U A V an d UGV trajectories: s min ≤ s k ≤ s max , (10) for all k ∈ J N − 1 K , wh e r e s min ∈ R > 0 and s max ∈ R > 0 denote the lower and u pper b ounds on the time duration between consecu ti ve time stam p s. Furthermo re, the UA V tr a je c tory and UGV traje c to ry a re coupled thr ough the U A V battery dynam ics. W e first consider the following bounds on the UA V battery le vel (measured in remaining flight time) e min ≤ e k ≤ e max , (11 ) X (km) Y (km ) r 0 r f Fig. 1: Exa mple p roblem with m A = 10 U A V task locations (red stars) and m G = 3 arms, who se associated UGV task locations are shown as magenta circles. for all k ∈ J N K , where e min ∈ R ≥ 0 and e max ∈ R > 0 denote the minimu m and ma x imum allo wed values fo r the battery lev el. In addition, the U A V battery dyn amics in c lude two modes: ch arging and discharging. The UA V can discharge its b attery at any locatio n, but it can recharge only when it lands on th e UGV , in which c a se th eir positions coin cide. W e model these require m ents using the follo wing constraints:      e k +1 = min( e k + κs k , e max ) , r A k = r G k , r A k +1 = r G k +1      ∨ { e k +1 = e k − s k } , (12) for all k ∈ J N − 1 K , where κ ∈ R > 0 denotes the U A V battery charging rate when the U A V is lan ded on the UGV . These constraints en sure that, between two consecutive poin ts along the U A V trajectory , the UA V ba ttery le vel either decreases, increases, or reache s capacity . In th e latter case, the U A V and UGV positions coincide at both trajecto ry poin ts. Moreover, matching positions at two consecutive trajector y poin ts imply that th e UA V and UGV position s coincide for all times between the correspo nding tim e stamps. In practice, we find that the constraint e k +1 = min( e k + κs k , e max ) is often overly restricti ve for n umer- ical solvers, pa rticularly when combined with the bou n ding constraints in (11). As a rem edy , we relax this constraint as e k +1 ≤ min ( e k + κs k , e max ) . (13) Intuitively , this relaxatio n allows the ba tter y level to be lower than the ideal charging outcom e, hence providing addition al flexibility for the solver without compromising the pr actical feasibility of the resulting solution. I I I . N O N L I N E A R S M O OT H I N G F O R D I S J U N C T I V E C O N S T R A I N T S The challen ge in solvin g th e trajectory optimization pro b - lem p roposed in Section II com es from the constraints. W e can d i vided the con straints discussed in Section II into two group s. The first grou p consists of constraints defin ed by smooth func tions: r A 1 = r G 1 = r 0 , r A N = r G N = r f , e 1 = e max ,   r A k +1 − r A k   2 ≤ v A max s k , k ∈ J N − 1 K , k p k +1 − p k k 1 ≤ v G max s k , k ∈ J N − 1 K , e min ≤ e k +1 ≤ e max , s min ≤ s k ≤ s max , k ∈ J N − 1 K , r G k = g ( p k ) , 0 m G ≤ p k ≤ p max , k ∈ J N K , p ⊤ k  1 m G 1 ⊤ m G − I m G  p k = 0 , k ∈ J N K . (14) The constraints in (14) are compatible with many algorithms for smooth non linear programs, such as in terior point meth- ods and augmen ted Lagrangian methods. The second grou p consists of disjun ctiv e con straints, given as follows: ∃ k ∈ J N K s.t. r A k = a i , i ∈ J m A K , ∃ k ∈ J N K s.t. b ⊤ j ( p k − p max ) = 0 , j ∈ J m G K ,      e k +1 ≤ e k + κs k r A k = r G k r A k +1 = r G k +1      ∨ { e k +1 = e k − s k } , k ∈ J N − 1 K . (15) Notice tha t h e re we changed the constraint in (13) to e k +1 ≤ e k + κs k . This change is lossless when we im p ose the constraints in (14) and (15) togethe r since the co nstraints in (14) already ensures that e k +1 ≤ e max for all k ∈ J N K . The disjunctive constrain ts in (15) pose unique challenges for op timization, as they r ely on logical OR operations and induce a disconne c te d feasible solution set. W e first discuss how to mo d el these con straints using d iscrete variables, leading to a mixed-integer nonlinear programming appro ach. W e then pro pose an alternative approach that first refo r- mulates the constrain ts in (15) as nonsmo oth constraints and subsequently a p proxim a tes them with smooth non linear function s tha t are compatible with nonlin ear p r ogramm in g. A. Dis junctive Constraints via Discr ete V aria b les A classical app r oach to model the disjunctive constraints in (15) using d iscrete binary variables. I n particular , we can reform u late these co nstraints as follows:   r A k − a i   ∞ ≤ µ (1 − U ki ) , k ∈ J N K , i ∈ J m A K , | b ⊤ j ( p k − p max ) | ≤ µ (1 − V kj ) , k ∈ J N K , j ∈ J m G K , | e A k +1 − e A k + s k | ≤ µ (1 − W k ) , k ∈ J N − 1 K , e A k +1 − e A k − κs k ≤ µW k , k ∈ J N − 1 K ,     r A k − r G k r A k +1 − r G k +1     ∞ ≤ µW k , k ∈ J N − 1 K , P N k =1 U ki = 1 , U ki ∈ { 0 , 1 } , k ∈ J N K , i ∈ J m A K , P N k =1 V kj = 1 , V kj ∈ { 0 , 1 } , k ∈ J N K , j ∈ J m G K , W k ∈ { 0 , 1 } , k ∈ J N − 1 K , (16) where µ ≫ 1 is a large positiv e scalar . The id ea is to first express each disjunctive constrain t as a finite set of candidate c onditions, o ne o f which mu st be satisfied. Next, we intr oduce binar y variables to encode the selection of th ese candidate conditions, and en sure feasibility by for cing the sum of the associated binary variables to equal one. W e f o rmulate the UA V -UGV trajectory op tim ization pro b - lem a mixed in teger non lin ear pr ogram (MINLP) . This progr am con tains the following variables n r A k , r G k , e k , p k , { U ki } m A i =1 , { V kj } m G j =1 o N k =1 ∪ { s k , W k } N − 1 k =1 . (17) W e formu late this MINLP as follows minimize V ariables in (17) P N − 1 k =1 s k subject to constraints in (14) and (1 6) . (18) W e can so lve MINLP above using branch- and-bo und–based methods combined with nonlinear programmin g algorithms. For details o n mod els, algorithm s, and prac tical solution methods for MINLP , we ref er interested readers to [18]. B. Dis junctive Constraints via Nonlinea r Smoothing One lim itatio n of the discrete- variable ap proach is that it leads to an expon ential growth in the number of p ossible values f or binary v ariables, which often makes scalable real- time solu tio ns impr actical. As an alternativ e, we introduce a continuo us mo deling app roach fo r the disjunctiv e constrain ts in (1 5) th at av oids the use of discrete v ariables. T o this end, we first refor m ulate th e constraints in (15) as follows: min k ∈ J N K   r A k − a i   2 = 0 , i ∈ J m A K , min k ∈ J N K | b ⊤ j ( p k − p max ) | = 0 , j ∈ J m G K , min   | e A k +1 − e A k + s k | ,       σ δ ( e A k +1 − e A k − κs k ) r A k − r G k r A k +1 − r G k +1       2   = 0 , k ∈ J N − 1 K , (19) where σ δ ( α ) =      0 , α ≤ 0 , 1 2 α 2 , 0 ≤ α ≤ δ, δ α − 1 2 δ 2 , α > δ (20) and δ ∈ R > 0 is a p arameter with small positive value. The key idea is to refo rmulate each disjuncti ve co nstraint u sing pointwise minimum of the violation of candidate condition s. Here, the func tio n ψ δ provides a smooth me a sure o f t he violation of inequality con straint. Note that the con straints in ( 19) are not compatible with algorithm s for smoo th no nlinear prog rams, du e to the n onsmoo th p ointwise minimum fun ction. A common approa c h to approximate th e poin twise m inimum f unction is via th e log-sum-exp function. Gi ven c 1 , c 2 , . . . , c n ∈ R , the appr oximation is a s f ollows min k ∈ J n K c k ≈ − 1 τ ln ( P n i =1 exp( − τ c k )) , (21) where τ ∈ R > 0 is a positive scaling parameter . The idea of th is appr o ximation is th at, as τ in creases, the exponential term corr e sponds to the sm a llest entr y will dominate the sum and cancel with the logar ithm function. Howe ver , th e log -sum-exp func tion often causes n umerical instabilities du e to the rapid growth of the e xpo nential func- tion. T o over co me this limitation, we propose to replace the log-sum- exp function a smoo th f unction constructed based on the ℓ p -norm . In pa r ticular , let c deno te the vector whose k -th en tr y is c k , we let c − 1 denote the elem entwise recipro cal of c , assuming c k > 0 for all k . W e can show tha t lim p →∞ ,p ∈ N 2 p √ n/   c − 1   2 p = 1 /   c − 1   ∞ = min k ∈ J n K c k . (22) Here the co efficient 2 p √ n nor malize the magnitude of the norm   c − 1   2 p . This limit implies that, for a suf ficiently large integer p , the ℓ p -norm p rovides a smoo th for mula to approx imate the pointwise m inimum functio n. T o av oid the singularity case where c contains a zero entry , we propose the following a pproxim ation f ormula min k ∈ J n K c k ≈  1 n P n i =1 ( c 2 k + ǫ 2 ) − p  − 1 2 p − ǫ, (23) where ǫ ∈ R > 0 is a small po siti ve parameter to ensure the function is well d efined, an d p ∈ N is a rea so nably large in teger ( in p ractice, p = 3 provide s a satisfying approx imation in simu lation). Notice that (2 3) reduce s to (22) if ǫ = 0 . Similar appro ximation ha s been used for log ical specifications in optim iza tion [1 9], [20]. W ith either o ne of the formulas above, we can app roximate the n onsmooth con straints in (19) as smo oth co nstraints and optimize the UA V -UGV trajecto ry optimization by solvin g a smooth nonlinear optimization pr oblem. This nonlinear pr ogram (NLP) con tains the following v ariables:  r A k , r G k , e k , p k  N k =1 ∪ { s k } N − 1 k =1 . (24) W e formu late this NLP as follows minimize V ariables in (24) P N − 1 k =1 s k subject to constraints in (14) and smoo th approx. of the constrain ts in (19) via (21) or (23). (25) Note that all of the variables in optimization (25) ar e continuo us ( i.e. , no discrete-valued v ariables) and all of the function s that appear in o ptimization (25) ar e differentiable. As a result, in principle one can solve this optimizatio n using standard alg orithms for no nlinear progra m s (NLP). Howe ver , the smooth nonlinear functions in (2 3) ar e o ften ill-condition ed. I n pr actice, we ob serve th at the augmented Lagrang ian method provide consistent ro bust p erforman ce for solving ill-conditio n ed NL P [21]. I V . N U M E R I C A L S I M U L AT I O N W e demonstra te the propo sed NLP mo del on a on e- U A V -one-UGV system with mu ltiple U A V and UGV task locations. Fig. 1 shows the route map. W ithin the NLP fram e - work, we com pare the perform ance of dif feren t smoothing function s and nonlinear prog ramming algorith ms. W e also compare the NLP model with the MINLP model to illu strate its scalability . A. Pr oblem Setup W e conside r a problem with m G = 3 , and th e UA V and UGV share the same in itial ( r 0 ) and final ( r f ) positions, as illustrated in Fig. 1. W e consider e max = 0 . 4 h , v A max = 36 km / h , v G max = 16 . 2 km / h , and κ = 1 . 5 . Additionally , we set e min = 0 , s min = 0 , and choose a suf ficiently large u pper bound s max = 10 h for the inequ ality constraints in ( 14). W e perfor m all simulation s on the Minnesota Superco mputing Institute cluster ( https://www.msi. umn.edu/ ), wh ich uses AMD EPYC 7702 pro cessors. Each simula tio n uses a single CPU core with 8 GB of allocated memory . B. C ompa rison of Differ en t Algorithms for No nlinear Pr o- grams W e compare the perf ormance of the propo sed NLP model in (25) using two co mmonly used algorithms for NLP: th e Augmente d Lagran gian Metho d (AL M) [21], [22] and the Interior-Point Method ( IPM) [23 ] . For the ALM, we use the L-BFGS me th od (imp lemented by fminunc in MATLA B ) to minimize the augmented Lagrang ia n (see [22]). For the IPM, we u se the fmincon fu nction in MATLAB with limited- memory Hessian appro ximation. N o te that the convergence of NLP algorithms is sensiti ve to initialization [24], [25]. This issue is particular ly relev ant in o ptimization (25), since approx imating nonsmooth functions with smooth one s leads to almost discontinuo us grad ients. T o m itig ate this is sue, we warm-start the NLP (ALM and IPM) b y initializing the UA V’ s and UGV’ s position trajecto r ies by setting their rendezvous locations as the p r ojections of the UA V task locations onto the UGV’ s road network. Fig. 2 shows the conver gence of the ALM and IPM under th e ℓ p -norm a pproxim ation in (23) and the log-sum- exp ap proxima tion in (21). W e compare the se methods using 100 problem instances with r andomly generated U A V task locations ( m A = 10 ). The constraint violatio n refers to the sum of the violations of all constraints in (14) and (1 5). W e set δ = 1 in (19), ǫ = 10 − 3 and p = 3 in (23), an d τ = 1 0 2 in (21). Am ong all the meth ods compared , the combination of ALM with the ℓ p -norm app roximation ach ie ves the best overall per f ormance. With the same smooth in g fun ctions, ALM consistently c o n verges to better loca l op tim a than IPM in t erms of bo th ob je c ti ve value and constraint violation . When combined with the ALM, the ℓ p -norm approxima- tion achieves lower constraint violation while m aintaining a compara b le objective value comp a r ed with the log-sum-exp approx imation. Fig. 3 sho ws the history of the U A V’ s re maining time of fligh t as well as the distance between the U A V and the UGV along an o ptimal trajector y (comp uted using the ALM with the ℓ p -norm ap p roximatio n) for the problem instance illustrated in Fig. 1. The shaded re gion s indicate periods during which the UA V is charging on the UGV . Fig. 3a shows that the pr o posed model in ( 25) allows both partial recharging of the U A V and covering multiple task locations during one discharging cycle of the U A V . ALM w/ ℓ p -norm ALM w/ log-sum-e xp IPM w/ ℓ p -norm IPM w/ log-sum-e xp 0 20 40 60 80 100 3 6 9 12 Time (s) Objecti ve Function V alue (a) Con verge nce of objectiv e func tion v alue. 0 20 40 60 80 100 10 − 6 10 − 4 10 − 2 10 0 10 2 Time (s) Constrai nt V iolatio n (b) Con verg ence of constraint violation. Fig. 2: Convergence o f different NLP algorithms and smo oth- ing functio ns o ver 100 prob lem instances with random ly generated U A V task locations with m A = 10 . The solid lines rep r esent the median values of the simulations, while the lo wer and upper bars ind icate the interquartile rang e, spanning from the 0.2 5 quantile to the 0.75 quantile. 0 1 2 3 0 0.2 0.4 Time (h) U A V Remai ning T oF (h) (a) History of U A V’ s remaining time-of-fl i ght (T oF). 0 1 2 3 0 3.0 6.0 Time (h) U A V -UGV Dista nce (km) (b) H i story of UA V -UGV distance ( km). Fig. 3 : U A V -UGV trajector y compu ted via NLP (ALM with ℓ p -norm ) fo r the p roblem shown in Fig. 1. Each r e d asterisk marks the point along the trajectory where the UA V reaches a task location . T ABLE I: Perform ance compa rison between the NLP model (25), solved using ALM with the ℓ p -norm approx - imation, and the MINLP model (18), with a maximum computatio n time of 6 hours. m A Method Succ. Rate Obj. Con. Vi o. Time (min) 2 NLP 100% 2.4351 5 . 0 e − 7 0.5 MINLP 98% 2.4352 5 . 7 e − 7 1.9 3 NLP 100% 2.4354 2 . 0 e − 6 0.5 MINLP 99% 2.4356 2 . 0 e − 6 18.8 4 NLP 100% 2.5188 1 . 9 e − 6 0.7 MINLP 87% 2.4352 1 . 9 e − 6 71.0 NLP (ALM w/ ℓ p -norm) MINLP 2 3 4 5 6 7 8 9 10 10 0 10 1 10 2 m A Time (min) Fig. 4: Median computation time versus m A over 100 Monte Carlo ru ns with r andomly genera ted aer ial task po sitions for each v alue of m A , with a maxim um comp utation time of 6 hours. Error bars indicate the 0.25–0. 7 5 quantiles, and the median constraint violation is on the or der of 10 − 6 . C. Comparison aga inst MINLP W e comp are the p erforma nce of the prop osed NL P mo del, solved u sing the ALM with the ℓ p -norm approx imation, against th e MINLP model. For the ALM, we use the same implementatio n as the one in Section IV -B. For M INLP algorithm s, we use the open-sour ce so lver Juniper [ 26] in the JuMP [27 ] fram ew ork. This solver im plements the Branch-an d-Bound algorithm with IPOPT [28] as the inner NLP solver (with limited-memo ry Hessian a pproxim ation). T able I and Fig. 4 compare the p erform a nce of NLP and MINLP over a varying number of U A V task locations m A , with a maximu m computation time of 6 hours. For each v alue o f m A , we comp are NLP and MINLP using 100 problem instances with r andomly generated U A V task locations. In T able I, we report the success rate and the median objective v alue, median c o nstraint violation, and median co mputationa l time. W e defin e the success rate based on whether the solver returns a valid solution within 6 h ours of com p utational time. In addition , we define the co nstraint violation for both the NLP and MINLP models as the sum of v iolations of all constrain ts in ( 14) and (15). Fig. 4 furth e r illustrates the scalab ility of the NLP model compare d with the MINLP mod e l. Overall, th e proposed NLP model ou tperfor m s the MINLP model as m A increases, especially in co mputation time. T able I shows that under similar le vel o f con straint violatio n, the NLP mod el consistently outperform s the MINL P model in terms of su c cess rate and computational time, with rela- ti vely minor increase in the objecti ve fu nction v alue when m A = 4 . Fig. 4 furth er shows that how the com putation time scales as m A increases (u nder similar le vel of constraint violation) . The p r oposed NLP model not o nly improves the computatio n time by up to two orders of magnitudes wh en m A ≤ 4 , b ut also en sure the computation time is within minutes when MINLP can not produc e solutio n s within hours as m A increases. These results showcase th e computatio nal benefits of the pro p osed m odel against MINLP . V . C O N C L U S I O N W e presented a non lin ear trajecto ry optimization mod e l fo r energy-shar in g UA V -UGV systems with multiple U A V and UGV task lo cations. W e based th is model on a smo othing approx imation o f disjunctiv e constraints, which elimin a te s the need for integer program ming. W e demonstrated the propo sed m odel on a on e-U A V -one-UGV system. Com- pared with m ixed-integer nonline a r prog ramming, this mo del reduces the co mputation time fr o m hours to minu tes in numerical simulation. In f uture work, we p la n to extend the current m odel to mu lti-U A V -multi-UGV systems and persistent monitorin g ap plications. R E F E R E N C E S [1] Y . Ding, B. Xin , and J. Chen, “ A re view of recen t advanc es in coor- dinati on between unmanned aeri al and ground vehicl es, ” Unmanned Syst. , v ol. 9, no. 02, pp. 97–117, 202 1. [2] R . Chai, Y . Guo, Z . Zuo, K. Chen, H. -S. Shin, and A. Tsourdos, “Co- operati ve motion planning and control for aeria l-ground auto nomous systems: Me thods and applicati ons, ” Pro g. Aerosp. Sci. , vol. 14 6, p. 101005, 20 24. [3] I. Munasing he, A. Perera, and R. C. Deo, “ A comprehen siv e revi ew of U A V -UGV collaborati on: Advanc ements and challe nges, ” J. Sens. Actuator Ne tw . , vol. 13, no. 6, p. 81, 2024. [4] N. Mathe w , S. L. Smith, and S . L . W aslander , “Multirobot rendezv ous plannin g for re chargi ng in persistent tasks, ” IEEE T rans. Robot. , vol. 31, no. 1, pp. 128–142, 2015. [5] P . Ma ini and P . B. Sujit, “On cooperatio n bet ween a fue l constra ined U A V and a refueli ng UGV for la rge scale mappi ng a pplicati ons, ” in Pr oc. IEEE Int. Conf . Unmanned Aircr . Syst. (ICU AS) , 2015, pp. 1370– 1377. [6] M. S. Mondal, S. Ramasamy , J. D. Humann, J. M. Dotterweic h, J.-P . F . Redding er , M. A. Childers, a nd P . A. Bhounsule, “Co operati ve multi - agent planning framew ork for fuel const rained U A V -UGV routing problem, ” J . Intell . Robot. Syst. , vol. 111, no. 1, p. 12, 2025. [7] P . Maini , K. Sundar , M. Singh, S. Rathinam, and P . B. Sujit, “Coop- erati ve aeria l–ground vehicle route planni ng with fuel constraint s for cov erage ap plicatio ns, ” IEEE T rans. Aero sp. Ele ctr on. Sy st. , vol . 55, no. 6, pp. 3016–302 8, 2019. [8] S. Ramasamy , M. S. Mondal, J. D. Humann, J. M. Dotterweic h, J.-P . F . Redding er , M. A. Childers, and P . A. Bhounsule, “Computati onally ef ficient m ulti-a gent opti mization framew ork for onli ne rout ing of U A V -UGV system, ” in Proc. IEE E Int. Conf. Autom. Sci. Eng. (CASE) , 2024, pp. 204–211. [9] T . Kim, A. P . V inod, and S. Di Cai rano, “De coupled trajectory plannin g for monitoring U A Vs and UGV carrier by reachable sets, ” in Pr oc. Amer . Cont r ol Conf. (ACC) , 2024, pp. 587–593. [10] Z . W ang and J.-B. Sheu, “V ehicle routing problem with drones, ” T ransp. Res. B : Methodol. , vol. 122, pp. 350–364, 2019. [11] S . G. Manyam, K. Sundar , and D. W . Casbeer , “Cooperat iv e routing for an ai r–ground v ehicle team–exac t algori thm, transformat ion method, and heuristic s, ” IEEE T rans. Autom. Sci. Eng. , vol. 17, no. 1, pp. 537– 547, 2019. [12] A. H. Eker , M. Oznigoly an, K. F . Karaagacl i, D. Gokalp, Y . Bickici , and F . Stroppa, “Coopera tiv e mission planni ng for multiple aerial and ground vehicles based on evol utionary computation, ” IE EE A ccess , 2025. [13] J . Diller , J. Sigler , and Q. Han, “Collabora tiv e path planning of energy- sharing drone-ugv teams for persisten t patrolli ng, ” J . Auton . T ransp. Syst. , 2026. [14] G. W u, M. Fa n, J . Shi, and Y . F eng, “Re inforcement le arning ba sed truck-a nd-drone coordinated deli very , ” IEE E T rans. A rtif . Intell. , vol. 4, no. 4, pp. 754–763, 2023. [15] M. S. Mon dal, S. Ramasamy , J. D. Humann, J. M. Dott erweich, J.- P . F . Redding er , M. A. Childers, and P . Bhounsul e, “ An attentio n-awa re deep reinforc ement lea rning framewor k for UA V -UGV collabora ti ve route planni ng, ” in Pr oc. IEEE/RSJ Int. Conf . Intell. Robots Syst. (IR OS) , 2024 , pp. 13 687–13 694. [16] X. Li, L. Y ao, M. Li, a nd B. Zhang, “Reinforc ement learni ng based colla borati ve path planning research for UA Vs and unmanne d ve hi- cles, ” in Proc. Int. Conf . Mac h. Learn. Intell. Comput. , 2025, pp. 595–603. [17] M. S. Mondal, S. Ramasamy , R. Rown ak, L. Russo, J . D. Humann, J. M. Dotterwe ich, J.-P . F . Reddinger , M. A. Chi lders, and P . A. Bhounsule , “Risk-aw are energy-c onstrained UA V -UGV coope rati ve routing using attent ion-guided rein forcement le arning, ” in Proc. IEEE Int. Conf . Robot. Autom. (ICRA ) , 2025, pp. 13 000–13 007. [18] I. E. Grossmann and Z. Krav anja, “Mi xed-int eger nonline ar program- ming: A surve y of algorithms and applicati ons, ” in Larg e-Scale Optim. Appl., P art II: Optimal Desi gn and Co ntr ol , 1997, pp. 73–100. [19] S . Uzun, P . Elango, P .-L. Garoche, and B. Acikmese, “Optimiza- tion with temporal and logical speci fications via generaliz ed m ean- based smooth robustn ess m easures, ” arXiv prep rint [math.OC] , 2024. [20] S . Uzun , B. Ac ¸ ıkmes ¸e, and S. Di Cairano, “Motion planning for in - formation acquisition via continuo us-time successiv e con ve xification, ” IEEE Co ntr ol Syst. Let t. , 2025. [21] L . Chen and A. Liao, “On the con vergenc e propertie s of a second -order augmente d lagrangian method for nonlinear progra mming problems with inequality co nstraints, ” J ournal of Optimi zation Theory and Applicat ions , vol. 187 , no. 1, pp. 248–26 5, 2020. [22] R. T . Rocka fellar , “ Augmented lagrange multiplie r functi ons and du- ality in noncon vex programming, ” SIAM J ournal on Cont r ol , v ol. 12, no. 2, pp. 268–285, 1974. [23] S . J. Wright, Pr imal-dual interio r-point methods . S IAM, 1997. [24] J . T . Betts, “Surve y of numerical methods for trajec tory optimization, ” J ournal of guidance , c ontr ol, and dynamics , vol. 21, no. 2, pp. 193– 207, 1998. [25] M. Y uan, R. J. Caverl y , and Y . Y u, “Filtering-l inearizat ion: A first- order method for noncon ve x traject ory optimizat ion with filter-ba sed warm-sta rting, ” in 2025 A merican Contr ol Confer ence (AC C) . IEEE, 2025, pp. 547–552. [26] O. Kr ¨ oge r , C. Cof frin, H. Hijazi , and H. Nagarajan, “Juniper: An open- source nonlinear bran ch-and-bound solve r in julia , ” in Inte grat ion of Const raint Pr ogrammin g, Artifici al Intellig ence, and Operatio ns Resear ch , W .-J. v an Hoe ve, Ed. Cham: Sp ringer Inte rnational Publishing , 2018, pp. 37 7–386. [27] M. Lubin, O. Do wson, J. D. Garcia, J. Huchette, B. Legat, and J. P . V ielma, “Jump 1.0: recent impro vement s to a model ing la nguage for mathemat ical opti mization, ” Mathe matical Pr ogr amming Compu- tation , v ol. 15 , no. 3, pp. 581–58 9, 2023. [28] A. W ¨ achter and L. T . Biegle r , “On the imple mentation of an interior - point filter line-sea rch algorithm for lar ge-scale nonlin ear program- ming, ” Mathematica l Prog ramming , vol. 106, no. 1, pp. 25–57, 200 6.