
CO TROL OF MOBILE ROBOT F,"J'J".H","'cJ....L
S1 G OPTI 12 TIO ETH 'D
By
KIRK ""ESSELOWSKI
Bachelor of Science
Kansas State niversity
Manhattan, Kan as
1990
Subnlitt d t.o th Graduat all g of
OklahOlna State niversity
in partial fulfillm nt of
the reqllir rnents for
the degr of
MASTER OF SCIE CE
Augllst, 2003
CONTROL OF MOBILE ROBOT FOR TIO
USI G OPTIMIZATIO METHODS
ThC'sis Ap\Howll:

=~_. 5_J_(3_/_~_O_3
Thesi" Advisor
/fti~t,T~__
~(
[)t',lIl of 1.111' (:raduall' (·"lIt,!!.<'
11
Acknowledgments
I lffer my t hanks for tllP guidaIl('t' or Dr. Rafael Fierro of th 1. RHE labora r
at Oklahoma State UiW'fsit,\',
111
Table of Contents
Chapter
1 Control Problems for MultiAgent Robotic Systems
1.1 Appli('(l t.iOllS fl)[" Svstellls of 10bile Robots
1.2 SWj)c 0(' Work fur the Th0.sis
1. 2.1 ()ft"lilll' Optilllizat ion for Coordinated antral
1.2.2 ~l()J<'! P1'0di('tiw COlltrol for Robot Formations
1.3 ChaptN SlIltlrtlarv
2 Logical Reasoning and Discrete Hybrid Automata
2.1 Li t.c'rals alJ( I Binary \ 'ari a hies
2.2 Dinctc' H.\'lJrid Alltolflata
:2 .2.1 Fill it (' Stal.<' \ [(lc!l iI)('S
2.2.:2 S\\'itdH'd Affilll' S.\'st.CltlS
:2.2.1 En'lIt ,c'lJ('ratlJr
2.2.l \Io<!c' Selc('{()J'
2.3 Chaptc'r SlIllllllar~'
3 Hybrid Systems in Mixed Logical Dynamical Form
3.1 COUV('X I1q~iolls ill State· Spit('(,
:3.2 \liXl'd Logical D~'lIaltli('al Fonll
.'3.2.1 AllalogtoDigit.al TIdatioIlships
:3.:2.:2 DigitaltoAIlaJop; TIc'liltiollships
:3.2.:3 L(Jgical l1C'iationships
:1.:2.1 ('olltillllOllS Dytl<ulIics
I\'
Page
1
2
G
I
12
13
13
IS
IS
IG
17
18
21
21
2'2
22
23
Chapter P
3.2.5 ixed Logical D:vnall1ical Form 2
3.3 Example 1: II Optimization Ov r a ll on pace 24
3.4 Examplr. 2: Iodeling a BOUllcillg Ball 26
3.5 Chapt(~r Summary 29
4 Solving Mixed Integer Programs 30
4.1 Th(~ Sirnplc'x Algori1hlll and Lin0ar Programming 31
4.2 BralldlamlBouud l\let.hocls for Example J 32
4.2.1 Optimization of a Linear F\Ine ion in a on orln/i x p cc 3
4.2.2 TIl<' SimpIPx Method Applied tu a SubProblem 3
4.2.3 Exploring the Optimizatioll Tn.>(· 37
4.3 using CPLEX to SolvE" Example '2 40
4.1: Chapter Summary 49
5 Control of a Team of Robots Using Optimization 52
3.1 The n()boFla~ Ganll' t':2
5.2 Rol)ot. 'lode'l Gt':
5.2.1 COllt.illIlOUS D.\'lIatlli(':,
5.:2.2 Fillitr St.at(~ l\Iachinc
5..3 The Rohot IVIodpl ill ifLD Form
·).3.1 A.llalogt.oDigital ReJa.tiollfihips
;).3.2 DigitaJtoAIIalog n(dat.i()I1sIJip;;
:'.3.3 Logical I1<'latio[lships
:) .3.1 D~'lI;Ull ics
~). 3.·) .\ 11 t (llllfl t.a
;).~~.;) COllstraiuts ou du' State 0(' tllC' D('f(~lId(']'
::>,4 GanH' Siul1llatioll n('sId ts
v
55
5G
57
GO
G1
G2
G3
G.3
Chapter Pa e
5.5 ChaptN SUlllmary
6 Model Predictive Control of Robot Formations 7n
G.l Pn'dicti\"(\ Coutrol ~lt'thods 70
0.2 j\/Ioclf'1 Pn'(lictiv(' oatrol Formulation t
6.3 \rpC [OJ' Cl Thn'('Rol>ot FormatioB (
G.4 SilIlulaJioll Results 0
6.5 Chapt('r SUflllllcU',\' 3
7 Dual1\1ode Model Predictive Control of Robot Formations
1.1 Lpa<!('rFollo\\'er Algorithms 84
1.'2 Probk'rn Forlllulatioll 8G
7.2.1 InputOutput Fe('u!Ja('k Lilil'a.riza ion 86
7.2.2 Dual Jod(' \Iou(·l Pn di('tiw Cuntrol 7
7..3 Stahility Analvsis 93
IA Sillllliatioll Re'sulfs 96
;.;:; Appl i('a t iOIl of MPC to Larg('I' FOl'1l1a t.i(JllS 98
7.6 Chaptl'l' 1J1ll11lar~' OD
8 Conclusions and Future Work 100
8.1 A F!<'xiblp FralrH'work for 1ultiRol>ot oordinatic)lI 101
.2 Futlln' \\'ork ill :\Iod('1 Pn'didiV<' Olltral. IgorithlIls ]().[
Bibliography lOG
YI
List of Figures
Figure
2.1 Illfonnatioll tim\' rrpresentation of a DRA.
2.2 A finitC' st.ate machilll' for a fourstat.e H. stem.
:3.1 A polyhC'dral partit.jolJ of state' space P.
3.2 A disjUllctioll of thrrshulds ddiuillg HlIonCOllvex region.
4.1 Optimizatiol1 tree for examp}l' 1.
4.2 OhjC'ctivp. fllll('tiOll value and hounds for rxamplf' 1.
~.3 Simulatioll of til(' MLD 1>OllIl<'ing ball mndl'l.
4.4 Pro('l~dllI'al flow cha.rt for gpuerating 1:iolutious.
G.l PI(\viug fidel for tIw g,U[lP, wi th ill itial posi tiOllS.
G.2 A fillit·(, statf' lIJachiJl(' for all attack<'!'.
S.3 A trial g;auH' wi tIl 3 d('f(lll(krs illld a.ttackpI'S.
5.1 I\: ilH'lIJat i(' itllmts fur tilt' trial gaUl(!.
5.0. tri;t1 gallH' with inputs milliIIlizf'c1.
f").(j I"':iw'mat.ir iupllt,s for thl' trial game with inputs millilllizf'd.
G.l G<'TI<'ral fOrln of all l\IPC 1'('SP01IS('.
G.2 A 3rol>ot formatioll.
G.3 Enllnt ion of tlH' 3robot forJllH t.ioll followi ng a traj<'etory.
GA Evolll t iOll of tIff' 3robot format.ion wi 1'.h no t.f'l'Ini nal cost appli('d.
6.,) Evolution of thl' 3robot formation frolll IHrg<~r initial l'ITorS.
/.1 Ol·hllitiow; for t hI' SBC cout rolll'I".
/.2 HpS]HHlS(' llsing the BC alld dllalUlod(' I\IPC a.lgorithms.
VIi
Pa e
12
1
19
19
38
39
;)0
S(
Si
GG
G7
68
G9
75
~
/
81
82
83
87
97

Figure Pa e
7.3 The dualmod algorithm with and", ithou coJli i u . ,t . 7
8.1 ooular ar ·hi e tun' for formation control. 10_
VlU
Nomenclature
CLF
CNF
DE:\.
EG
FS~I
I\UP
ILD
\1S
IPC
SAS
BC
Control L:vapulIO'· Functiull
COIljull<:tiV<' 0:01'l1lill Form
Di::;crd.(' Hybrid .\utulllata
EV('lIt G<'lIf'l'atur
Fillitc'Statr I\lachiue
~Iix<'dIlltcg('r PrograIll
i\Iix<'u Lugical D,Yualllical
l\loc}r Select.or
Iodl'l Pr<'clictive out.rol
Switch<'d Affine' Sy::;teIll
S<'para iOII Bearing COlltW!
IX
List of Symbols
HI) Pn'diclioll hurizoll
lL Li teral stat.C'lHplit
U Vc'ctor of ystem ill]Jut varialJk'
UIJ Vert.or of bill<1.l',\'yaluc'd illpUt. variablc's
u, \.'('dor of l'I'<1.Ivaltl('d input nl.riablps
U Scrips of ill plI t varial>1(~s
U" Optimal seric's of input ,'mial>l('s
U Fl'asil)lp s('ries of input YC:lriables
Uo Optimal feasiIJ\c' seri('s of input Huial>les
i:iU Series of iu('['('ulC'lltaJ input variables
U u]Jsc·t oj' JR." ill "'!lith SOllie U Illust reside
l' Kim'lIlatic ill)Jut liuear V(·locity for a I'Ubot
i' Ob,iecti \'(. fUllctioll
x \'('c'tor of slat.e "ariables
XIJ \'ector oj' hillClry\'allleJ (I3oolc'all) st.ate variablps
X 7' "('dol' lJf l'C'al\'ct!lled (collLiuuous) state' variaLles
X Sul>s<'t, of JR.71 ill which tiOllJ(' x BlIlst reside'
y \'('dol' of systc'lll out.put \'ariai>lcs
z \ '\:,cl or of auxiliar:v n~al\'alupJ variablc's
t5 \"(~('tor of <1.11xiliar.\· Bookall variable'S
Ii IucrcllH'utal \'a.luc' of ouj(,eti \'(. hlllctiou
.....: KiIH.'lllatir illpUI augular \clarity for a 1'OI>ot
x

Chapter 1
Control Problems for MultiAgent
Robotic Systems
III this chapter we' defillP what we hope to achirve by dr\ doping optilllizationIJa..<;c·u
control nwthods for multiagent robotic .~ st~ms, and motivate our approa'h b.\' discussing
past efforts presented in the 'outro] literature. We introduce thp two problem~
for control of teams of autonomous rnobih' robots to he iuveHtigat.r.d ill tlli~ t.Iw~iti:
a cooprrativr rontrol problem solvrd Il~ing oH:'lim~ optiJllizatioll. and a [oml Lt.io[l
control probll"m solvC'd using recur'iv(!, recedinghorizon onlim' (but no(. I'('altiHw)
optimization. Both types of syste"Ills are sirnulat.ed aud discussed, but llO verification
using real rohots was undertak('n.
The C'ooppl'ative control problelll is solved by modclling H team of Illouile rohots
as a h~!hrid svst('JlI and controlling the systmn by solving it rnixf'elint(·gm program
offline to find some suitable inputs for the system. Th~ formation cont.rol prohlc'lIl is
solved using n~(,f'dingllOrizonmodd plwli('tive ('outrol 1Il(lthods, and \VI' COlllpan' aIlel
cont.rast these methods to somf' previolls work Ilsing passive' state' f(·(·dback coutrol
laws. \Ve rl iscllss some difficu Itics ill optimizatioll, amI poillt 011 t SOll1(' potl'l1t.ial
dra\vbacks to model pI'f'dictive ('outro] that w(' will 11lwl to address ill the> cl{'wJOPIlJeut.
1

\Ve justify the inclusion f h e two di,}. ax ar Ll1
how concep, from ,ach may e ,ntuall. b ;pplie lad, el ping a fram ;w rk £ r
online recedillghorizon 'ontr \ of eam. f au 110m 11 I' b t. m d 1.d h. brid
stem', These idea: are intra luced in thi. chapt I' and further e panel ,d , h 11 P sible
future work i di, eu ,ed ill hapt I' , H wevrI'. rl'gardless of an, future work
ill this area, the main concept, introduced in thi. tli si ll, iog optimization me hods
to coutrol h,vbrid . ystem . and the potential ben fits and hallengcs associated
with recedinghorizon. optimization ba,srd GOlltrol m ,thodsmay have appli ati us
in man:v exciting areas 01' control that are Hllrler <if' elopment. (spe, for example. [1]
for a survey of current trrnels),
1.1 Applications for Systems of Mobile Robots
There are applications wherp a coordillatcd tcaw of lIlohile robots is mon' suitable
than a siugle robot, especially whcn~ a. distributed system of sl'nsors i' advantageou '.
For t'xamph\ iII sc~ar('balldl'f'st'ue opera tiOllS. l.!('ploYlI H'I 11 of IWUl r robot s llV('r all
arca cau allow fur Ilion' thorough ami faster cuv(~rag{'[2], Ot \1('1' apl'l icatiolls an' in
geographical cXlJloratioll. lIlollitorill~\ and mill(' sw('cpillg [3]. Yet without'. nJonlillatiug
the IlHlV('IlH'llt of tiLl' ageuts, any advilllt.age of Jl1lllt.irubot d('ploYlIl<'lIt'. lila I>e
lost auel dalllagillg collisions or iutl'l'fereJ)(;(' lllay (wcur.
ContrullJrobl('lOs fur Ulultiple autouOIllU\lS vchides call haV(' aPIJ!icatioIls ill a wiel('
\'ariet~· of fields. At, OIl(' ('lid of tIll' SIH'l'1l'1lUl are 1II1lltiall;cllt S,YSt,{'lllS f bat art' UOt.
working together. Imt Iwed to 1)(' coordillat.cd to tJH' extcllt that tlH'y do not. iutpl'fcr('
with our allotheL Typical Clpplicat.iolls ill t.his area an' iu airtraffic tollt.1'01. At 1.11('
oth('r Plld of the SlJl'ctl'llJll arc tcalllS of rohots olwl'atiug roo(><'l'ativdy alld 1ll0\"iug ill
forlllatiolls wit It pn>cisely defiued g('ollirt ric's that IIc('d to rl'act to thl' C'll\'iJ'OIIllIt'nt
to <\Chic\'(' st.at.ed goals. III the middIc' gr() lIuli tIH'I'(' an' gru\I ps of \'('h ides that have
2

1'1 n
h VI r ha
me appli and
on ider fi lU hi. hI' d I' p. rum
m.
orne common goal , bu m, n h
ontrol problem in this ar a i
or herding b havi I of anim 1 .. \
. how ho\\ the con 1'01 problem
f multiple autonomous velli 1
III conflict r' olution prabl m., traj ct Ii < are planned for vehicles operat.ing
independently in the same 'Pace but with poten iall~ onfiieting goal, ircTaft. trajectory
planning ha.·, been, tudied b, [4], [5] a.nd [G] using mixf'dintrgcr lillrar programming.
The work in [7] us s th ,d Darni window approa h to a oid obstad s, and
incluc!c=.s the stability analy. is. The C:' methods often a.dopt w r, 'implifi.d model. f
\'ehidc behavior, but pres nt int"resting application for mixedintep;er progntrnmin
fIle tllOds.
In cooper'ative control 'PT"Oblem' vehicles mov in a coordinated fashion to achir
some commuu goal and/or seek to maiutain some grometrical relationships alllollg
themselves. Often mOVp.rnellt is dictated by measurement of gradients of some actual
sensor 1IIeaSnrelIll'nts, or some art ificial potential field, 01 11 tioHS dofill<'d wi tit illterrobot
distance relatiom:;hip::; are f'xplored in [8], wiIprp mc Ito Is to Ilwasu['(' alld projcct,
gradirllt illforrllafion are discu sp.d. The. applicatiolls for th .Sl' mf'tborls arr in. for
exarupll', dat.a acquisition in liUgl~ an'as such a ocealls where' the 1l10St. advantap;eou
arranf!/'1IJ011t of sensors lJlay not be to distribut(' t.hem eV(,llly, hll t t.o haw them
adapt to COllcentrate more sensors in areas wherf' the measurod v riable has steeper
?;raelients.
An optimil.mtion approach is used for path plalllling auel collisiotl aVl)idance of
mariIll' vehidr.s in [9]. For a vehiclr with multiple desir/:'rt bl~haviO['s, such ali collisioll
avoidance and path lellgth Illillimizat.ioll. UtI' aut.ltors <!('nlH' a wult ipkobjl,('t.iv(' fUllction.
wi t lJ a different itllalytical expn'ssion for each hrha.vior. Tl\()~1 develop .'It rat.pgil's
to pan' dOWll the' sparch space. and solw an intf'rval pwgramrninp; problpIll rrcllI'sivdy.
imilar appli iOQ i lJl ,b f l' b ti gam uch . Rob UP.
occer. .L\1though h imm dia . appli aLi n i. impl~ an in rc lIe ia t ~ r
robot team' playing occer, he competiti II prO' ide' a platform for h0 dey I pmellt
of robotic control, ensor, communi 'ation and mechani I yst .ms. ue 'e sfnl
('antral strat.egies are often Lehaviorb' sed with robots operating in certain mode
"vith hE.uristic rule's determining how cadi robot f(~acts to best lWIIcfit it:i 1",('<\111. \"rt
withill this class, good control 8trat(~gies rna, require morr thau impl, a (,t of "ifthen
else" rules defining th mode wit hing.
Coordillated for'mation control refer to trajm:tor r tracking h a formation of
robots having a particular de ired geometry. Thi is has application' where rigorous
coverag(, of an area is required, snch a'i ill grographical exploration or min
SW<:'r.pillg[3]. Conceptually. fighter jets Hying in a delta formation at an air how are
an example of (rnanuaIJy coutrollcd) autonomous vehick' tracking a trajertor.v in
formation,
ThE'rr have beE n a variet.v of approaches to formation control; wr can makE' one
distillct ion hased Oll tht' (lPgree of CE Illl'alizatioll, The 1('ath~rf()lIower al'prmu:lws
an' ('('lltraliz(·ld LH'causf' they involve a t('anl of rohots followiug j\ r('al INHl<or robot
or j\ virtual robot or [orruatioll of robots. OlJ(~ approadl for ('oul.wl of fOrJ/lat.iolls
of robots is ru USf' all input out]J'lLt feed/mcA; Liw~u1'izuJi(m COli t wI law appli('d to it
"follower" rohot. at SCHU(' spccifif'c! wlat.iollship tu a "lcad!'r" rollot. [IO] [11] . 1 siug
this t.(I('hniqu(\ cOllJplex fonllatiolls cal! he achieved \Vitil ddiued n'!aLiolls for chaius
of 1'0 hot.s. EadI robot ill th<.' format iOIl III ay 1)(1 followi Il I!: onf' robot alld I(·;\( Iilll!: OU('
or lllon' othel' robots. Stabi lity of forruatio/l ('ontrol l1si nil; t.his lIIf't,hod ('all 1)(' proVC'1l
lIwkr appropriate conditions [12], 0111<'1' n'srarchrrs ItaV<' IH'O"idl'r! SOllJ('what sit/lilar
algorit hillS llsillg \"inual lead(~rs [13].
Lpctd('rfollowPl' approadlf's ma." 1H' SIlS('(')H iIlk to faihJr(' of the kad (II' or l('ad('l's.
or ma.,," sllff('r from a lack o{"fC'pclba('k to till' kadn if t.11(' followers aI'(' i/lljH'dC'd. OIl tlw
other haud. he man' de .n raliz. d p imizationb , . m h d b morc mbu t
aut! man' fil'xiblf' ill rnC'Pting p rforrUalH'(' rrqllireruf'ut. . bllt may 1)(' COlllI utatiullall '
expclIsin'. A d(,(,(,lltralir,pd approach is takru in [11] ",hrH' the' lack of robu,tll('s: of
hr l<~ad(')'f()llow(~rllI(·thods is ov('r('01l1(' bv using dev '1' elclllf'lltary formatioll lIIoves,
Fonnati()u ])('havior is drinm ill a '('ric's of fonllatioll moves that an' assullH'd t.o be
fcasiIJlt' aud ['('snIt. ill no collisio11s.
TIll' ("uordinat.iou of upperIrwl or SHIH.'n'lSOr,V control with low0rI(,\,('1 or local
control has hl'('ll all important topic for formatio11 cont.rol [IG] [16], III [15], (,(llltrol
graphs an' ddiIH'd t.hat allmx;
the scparat.iOII of cout.rol task.' int.o group fonnat.ioll
auel local ("ol1trol I(~\'f'ls, A h('lll'ist.ic algorithm is propOS('<\ iu lieu of a, mixC'dillt.(·ger
program t.o to opt.imiz(' the fOI"matioll ('Ollt 1'01 1('\'1'1. III [1 G], st.aIJili t.y is showu for
control of fulh'act.ni:1tru (lJoloIlOlpi(') n~hicl(~s, usiug a lll('t.ltod wit.h "irtllal I('aelers
and potrlltials, Howc'vpr I thrs(~ llH'thods cI () lIot sh(m' a fkxi1,Ir fraI1l('\o\'()l'k for (,olltrol
at tllr higher lrn'l. It. way OP possible to devdoJ> ulOdrls of t('ams of rolJOts as
iJ,vbrid systems. and to implement cOl1trol scIlf'lllP:'j usillg l'('(wlillghorizoJl pn'di('tiw
cont.rol llirthods. SIl('('C'SSflll C\PV<'IOplll('lIt of thl'sl' llJ(,tlJOds lila," c'('rtailll,v d('I)('l1d 011
au lllld('rstalldiug of all the' areas illt.r()<!U('('<! iJl this t.he'sis (alld IIHIIlY III ()['(,) , III tllat
SE'I1S(', future' work I1ta~· tie tugetl!n tl)(· two tyP('s of silllldatiolis !H'('S('II('d ill tllis
thesis,
1.2 Scope of Work for the Thesis
hi this work we ",ill addl'(,ss two Illltitiag('lIt ro],oti(' s:vsU'J11 c'ollt.rol pro],lc'llIs (hat. fi(.
within tIll' spectrulJl of (,olltrol [lW"IC'I1IS c!l'snilH'd al)()vc~,
;)
The first problem we will d moust.rate i a co rdin t d c n r 1 r bl ill v hen' WE"
develop a control. y. tern fm a • iroulation of a rob i game PI' po .d by [17] call ~d
RohoFlag. Most uccessful control ystE'll1 for th R b up rob ti liO er ha (> been
simple and reaetiVf'[18], Howrver, the rul s of RoboFlag 801' mol''' omplicated, so t!l
strategic planning should be rewRIded. Thi i intended 0 rno iv t developmt'llt of
control systems Oil a supen isor:v l~'v<~l.
As Ro!>oFlag is quitr a complicat d, we' pre. ent a 'imulatioll of ouly one aspect
of th(~ gallJe, ill the manner of [19], \1\ f' model the team a,"i a hyb1'id sy ·f; m, that
i . Ollt' having both eventdriven and tim driven d~rnami(:s. \lVe sol.ve tlH' control
problem by formulating an optimization problem and solving it off liIle u'inp; mixed
illtegpr linear programming, EVE'1l if HO RoboFlag ('orn!wti (iOllS ('V(')" dp\,r!op, t!lis
work df'fllonstl'ates how we can develop a rigorous framework for suppl'visory control
of morf' general multiag nt robotic s~rstcms. Eventually, iIllpkrnelltatiuli of oilline
rnodpl prpdict iw ('ontrol for hybrid sy trms may stem from such a frampwork, possihly
luwin!l; applications in a wide variE'ty uf an~as.
\\"(' dC'vrlop Ol\I' descriphon and anal,vsis of the game hy first iutrudnrillp: logical
r('latiollships in Clmptf'l' 2 and describing Iiow we' can repr seut tlwm l1sinp: S(,ts of'
simultaneuus pquations with discn te 1>inar,V variables, One way to d(~s('rihr a hybrid
syst.em is h~ IIsillg a disC7'ete hybr'id automata (DHA) we' ddinl' the DHA and sbow
how it nUl 1)(' represented ill rniJ:ed loqical dyrwmiml (MLD) form in hapter ~3, [(
oH'lillf' simulation of a DH in 1\ILD form results ill a miJ:(~d inteyer p1'Oymm( MIP)
opt,ill1izatioIl problpm. We show bow" UP's can be solved in Cha.pter 4, and !JH'1l in
Chapl(']' .j \\'(' use t.h\' RohoFlag gamr to dprnollstrat(' a lllO[,(' illtcrt>stiIlg ('xall1plc of
tlwsp IHf't.hods, As a. product of this, we have contrihutl'd software> that ha."i ..dlowc<!
11S t.o usp t 1)(' appropriate commc'ITial software packages to salVI' similar problelIls.

1.2.1 Offline Optimization for Ooordinat d Con rol
G
1.2.2 Model Predictive Control for Robot Form tion
Our secoILd problem is a traje('tor~' tracking format.ioll cont.rol proLlr.lll, when' we
show ('Ollt wI s~'stellls for trajec tory t rackiug of nllli t.iagE'llt. robot format iOIlS, III COlltrast.
to t l)(' offlim' silllulatioll IlsC'd for the nohoFlag gamC' our sillllllal'iolls ill this
arca apply n'cursiw olllilH' 1/l.()df'lpn~didive control ( IP ) algorithms to cOlltrol
[ormat.iolls of robots, Geuerall~', 'fPC algorithms rel~' Oil an optirnizatiOlI of a predicted
mockl rpspOILsr wit.h rpsped to plant inpllt to uPt,cl,'min<' the brst illput dlang'l's
for a giV<'1l statp. Thl' ~IPC algorithms call ill geJlrral handle 1I0nlillrar Illo(kls, alld
lllay 1)(' l'('COllfiglllabk, Eithpr hard (,ollst.raillt.s (that call1lot bC' violatl'd) or Sort. C'OIIstraints
(tlta t. nUl 1)(' yiolated 1,u t ,,'i tit .'lOllI(' pl'nal t,\") call he illcorponHPd into the
optimizatioll, gi"iug ~IPC a pou'n1.ia I advallt.ag(1 over passive statr [('ed hack ('0111.1'01
laws ill SUU!P applicatiolls,
RI'c('utl,v, Dunbar and Murray [20] Iti"l.Vl' ckvdo))()d a '1PC control algorit,lun for
tra.j('ctOI'~r tra.ckill~ of teams of robot.s in certain forIllations, This al~orithw lISI~S a
t('[IniIlnl cost. to I'llSlJI'(' sta hi! ity, all< I rp}jl's OIl s,n 11Ill('t.r~' of a tri;lll p;lIlar ron Ililt.ioll
to rq~lllatl' iut.erroho{ sl'paratioI! ..we! distaw'(1 "WIll t.1ll' t.rajl'ct.ory 1'<'1'1'1'('1[('(' POilll,
Othl')' [('C('ut work ('hat IIS{'S t,erIlIillal costs to ('IISIIl'C' stability for 'l1I1'OllUII10llS \'('Ilicks
lIa,s 1)('('11 l:outril)]Jl'('d b,\' [21] all<1 [22], III Chaptl'!' GW(' PJ'(~S('lIt il SillllllilLiolJ (If a tltn'I'rohot
fOrtluHiOIl l'ollt.l'Ol prohklll IlSillg an 1PC algorit.luJl ill t.hl' mallJlI'1" or [20], and
disCllss what COHdit.iO[lS lllay ))(' IJ('cessary for st.abilit,v, \V(' do HOt. COIlsid('1' ;1I1~' [ol'lllal
proof of st.ahilit~': lllstl'ad WI' discuss what St.I~PS iUI' 11l~('('SS;UV ill ordl'r t.o l'IISIIl'l' llial
thl' sillllliatioll I'('SpOHSI' looks till' ",a.\' WI' illl('lId it t.o I)('haw, \VI' d('U10lISt ra 1.« I that.
l'('lllO\'al of the 1<'nuinal ('ost does J'(~sllit iu a S~'St<'IJI that has llIlstal}I(, d,"U<lltlics,
As rlll altel'natin' to lIlC'tl!oe!s lIs'lUg tl'J'lllillrd cost.s, ill Chapt.C'l' 7 \\'(' jll'l's('111 a
TlwtllOd llsing j('rmillal (,ollst.raiuts t.o 1'1lSIII'l' stahilit,V or thl' i\1PC algoritllllL, From
th(' tel'llliw,d coustraillt rl'l:?;iou, \\'(' I lSI' a sl abilizillg local ('(Hltrolh'r ill th(' Illrlllll('r of
[23] aud [21], lI'sllll illg iu a duallll,ode ~IPC algoritlulJ, This llI(\~' hm'l' thl' adnUlt.ag('
7

of concept.ual ·implicit. but rna, uffer th eli ad antag 1. ptimiza i n
prablem~ outside a Bec .. saril Limit.d r ri n of c n crgence £ rhe al ri hm. lnd .d
we of1'01' a proof of stabilit, for thi alg ritlul1. although it requir . umptioll. hat
may bE' difficllh to meet in practicC'. \t\f contI" t the dualIllode' MP coutrollrr \ ith
tlH:' use of an inpu output feedback linearizat.ion controlla'w alolle. In the dualmode
iPC algorithm, we cau im:Ol'porate a colli. ion co t t€'rm to reduc€' or dfrcti f'l.
elimiuate the possibility of intf'rrabot collision.. Wf' disrllSs the ad antages all I
eli a.dvantage of dualmod 'ontroJ versu. 0 h r pr die i e ontr ll~r , and diSCUSA
applications for ('ontral graphs of lar er arbitrary formations. The concept. allC l
anal.\ is in Chapter 7 an. derived fr rn ideas in [23] aud [24], IlIlt we are uot a"varr. of
t'llE'ir PWViUllS applir.atioll to this particular problem[25].
1.3 Chapter Summary
III t his chapter we haY<' illt.roduced the control systems we will analyze and simulat ,
aud :,>1Iown how they fit into the broader Hpe trurn of Illultipll' autonomolls vehick
('outwl problems. We will introduce a.nd 'olv(' wo broad t..vp. f IlIulLiple robllt
control prohlplns: the RoboFlag ganle i' an example of a coop~ rati e coutrol problem,
and formation trajex:tory tracking if; an example of a formation . ntrol prahl fll. h
of t.hem are simulated using optimizatiouba.':le I methods; thC' former i. Holved in all
offline simulation, and tbE' latter ill a recursive manner that would f'wntually be
siIlllllat~donline IIsing model predictive ('outrol met.hods. By dewloping cOlltrol
tra~ep;ips ill thesf' two area",;, w(' hop(' C'vt'utually t.o dpv<'!o!J an lIu<!cl'stalld illg' of
('Oll<'l'pts that. will allow lIS to CTPatp a fkxiblc' awl ('(Jlupn'h('lIsiY<' frallH'work for
C'olltrol of rohot teams 011 thr supervisory Ie'vel and thl' loral kvd.
8

Chapter 2
Logical Reasoning and Discrete
Hybrid Automata
The supervi 'or. level of control can be 1'('prf'sent d a:::; a set of logical ruks thnt cue'
analogolls to all intelligent operator d(,tining tilr tasks of tlt(· ImV<'r leV(~l d.vmutlic
control::!. Qne example of a logical rult' is the famil iar "if then' statrlll('ut, bll t w('
show that all~' proposi tioBa,) logie statl'IIl('lI I can 1)(' ['(IIH'C'SC'IlI I,d ct.<; a ['(·Ia tionsl.i p 1)('twtJpn
biIlary variablC's, discllssiolJ of n()oh~all logic aBel lillear t.hn'sholcl cOllditiollS
prO\'iops till' COlltl'xt for a hrief. informal introductioll 1.0 It~'brid s,vst,(IUlS. TIJ(I ['ortlla!'.
that WI' will uSt' to describe hyorid s,vstC'nlS is the distTet(' h~'I)['id automatoll, which
wr will ddinr. ill term. of its major components, sprcificall,v t.lll' finit state machine,
swit.chpd affine s,vstem. mane selector. all d ('Ve'llt g('lIprato[', All ulld(~rstand iug of
tllf'sf' eOlllpOlwllts will sd, the stag(' for a syst.ematic Illptltod ['ur 11IOdr.lillj.!; II,vhrid
s~'ste)lJs ,
2.1 Literals and Binary Variables
A litfml is any logical statemrnt that is citltrr trll<' or falSI'. It is COllv(~llient to
a,"'isociatl' some' binar,v \'(uiabl!'. say fJ E {O. 1}. with som(' lit(']"al L. so that () = 1 +7
9

L = fTlL(' alld 6 = 0 H L = fal e. , e rna, al 0 ita C:' a literal If'fin d in erm
of ot hpr Ii terals as a 80 lean expr s. ion. ne commou \ a. t repre C:'llt
expression is in conjuncti1 e normal form (C F), whi II cau be l:'xpm . t>d a.
olean
where 1\ represent. a logical "and" operator and V )'('pres llts a logical '01''' operat. 1'.
Ea,(']1 of t.he expressious (.1 11 V 0I 2 V ... ) is rallf'd a disjunction aud the illt.ersection f
two Boolean expI'0.ssions (.) 1\ (.) i called a conjunction, so th) N i' a couj unction of
disjunctions. AllY Boolean expressi n ran bE' cast iuto Cl F by using the appropria
rela tions[26].
V\'p fall l1lodel logi('al (\xpreSSiolls llSillg the associated binary variable, by ('011
stl'Uet.iug' the appropriate serif's of expres, iOlls. The simplest expressiolls for litNals
(2.2)
and
(2.3)
with6i E {O.l},i = {1.2,3}. hlordt'rloilllpklll(llltalin(~arprognulHlIillgS()llltion
to logical problems, we wOllld like' to cast these rclatiOllsllips into furms suggested by
[27], so that (2.1) is eqllivalpnt to
6J + O2  e5:1 < 1
6) + 0:1 < 0
e52 + e5:1 < () (2.4)
alld thp logical "or" of (2.3) is equivalf' III to
10

6]  02 + 03 < 0
01  O:i < 0
/)2  6:1 < 0 (2.5)
Similarly, it. ea. y to spr thr logical "ifthell relationship Ita. h qui va!t· 1)('('
whefr I is rr(1<1 '"not L" auo illoicates t.he logical iuvrrse of L. Thus, L f1'I/(' +t
L = fal.,;('. This ('xpr<>ssioll is simply <,quivaknt to 6,  62 ~ O.
("Olllplrtp till' list of rxpn'ssiollS for proposit.ions as follows [27]:
imilarl.v, W(J nUl
L, 7 £.2 ¢:} 61 6:l < 0
L 1 +t £2 ¢:} 6, (h 0
£, V L'l ¢:} 6, + 152 > 1
L, 1\ L2 ¢:} 61 1, b'2 1
£1 ¢:} 151 0
LI EEl L2 ¢:} 61 + ()'2 (2.7)
sill/!; the ('xpn'ssions ill (2.7). we ("all build an.v logical r('laJiollsliip. For exalllpl(',
ill (2.G) to o!JtClill
15 1 + 62 + 6:l <
61  6:\ < [)
6:l + (h > (:2.8 )
III a similar IllauII ('!', \\'(' nlu (!<'vrlo)l ilH'qllaliti('s rxpn'ssillg all)' dpsin'd logi('al n,lationsiJip.
For a full truth t.ahle of rdatiotlslii[>s. s('c' [2G]. L"sillg t hc's(' ]'('latiollSliips, \v('
('an sol \"(' logind ]'('asolliug prol>lc'lJJs IJ.\' <!('\"('Io[>i Ill!; ill teg('I" programs to sol V(' ('qua t iOlls
listed ill tll<' f01"l1l or (2.4) and (2.,)), as showll iu [27].
11

Event Generator
(AID Converter)
Finite State Machine
Switched Affine System
Mode Selector
(D/A Converter)
Figlln~ 2.1: A symbolic illfonlliltioll How rcprcsf'lltatioll of a discrete hybrid automatoll.
2.2 Discrete Hybrid Automata
A hylnid system is Ollc' that. has both (,\,(,lltdri V('II auel tilllPdri vcu ([."uam ics. We'
will dl'scril)(' a h~;brid system I>~: llsing a mixedilltcger state v('('tor x = [Xl'I' X""']'1 '
\\' h('n' x,. E JR." is a v('('t,Ol' of 1'('(\1 vallied statl' variablrs and XIJ E {O, 1r" is a \'('('(,01'
of bimll"\·valul·d stall' varialJks. For Im'vity, WI' will n>('('1' til til<' ('OIIlPOI)('II1.s of Xl'
as real \·ariahl(·s Ol' [('al statp variahl('s, ami tiH' l'Oll1pOIH'Uts of Xb as biHar:\' \'ariabl<~s
01' billar.\' stat(' \'ariabll's !lrI'l'iuaft(,I', OU(' ~'ay of [('pJ'('s('lJtiug a hybrid s.\'SI(·1I1 is as
a rli.w;rde hybrid llut07fwton (DHA), whit'll consists I,f' 11 fini/;e ,c;tlLt(~ nw,r:h.irw hcwiug
biuary stah' variabl('s awl CI .'iwitched a,ffinr~ .'iY8tem having I'('al staU' varial>ll's [28].
Thf' fillit(' statl' mac:hiw' and swit,dH'd aftill<' s.vst.c'tu an' J(·lat.c'd I>y all event f!r'7U'1n!o'/'
alld a 71/,adr: selector' as shmYJI ill Figur<' 2.1.
TllP DHA fOrIlllllations arc' so]n'<! ill disl·n't(· tillH'. bill lIJ tlI<' 1it,('ratllI'(' or 1J~'hrid
\vst('lllS iIlt('w'rvalw'd sl'at(' variahl<'s ,U(' (Jfl,(~l1 call<'d "dis<T('l,('" va ria bl(~s. :\ dis('J'('t,(·
varia!>\<> cO\Lld take illtpger ndlH's alld d(ll's llot h(lv(' to h(' I>illar~:. but ill t.his work \v('
Oltl.\ ("oIlsid('r l>iIlalTvallH~d dis('1'('tl' variahlrs. To avoid ('(JIlfllSioIl. ill this ['C'POl't WE'
12
will alway u e the word' binar. ' to rrpr .. elit autom' ton e and biua! . riable
and "dis rete' to rE'present discret.r time. imilarly ill tlu literatUl'r often th . real
variables are called "continuous' ,·ariable. ,
2.2.1 Finite State Machines
A finite state machinE' (or finite' alltomaton) consists of a finite 'tate' spacr sd aml
a transition rule Se:'t whosr ('killC'IltS are calle:'c1 edge.c;, 1 he:' edges dC'fille ('vpnts that
dictate transitiollS uetw('u} statps of the' fillil e statr Illachin . The eclge. definp tlJe'
state of the finite state machinf\ accordillg to SOllle Boolean rule:'[28]
x/l(k + 1)
Y"l.k)
fU(Xb(k), udk) e5{.g(k))
g13(XIJ(k), ub(k), e5"g(k)) (2.9)
wh(~re xdk) is a Booleau state' variable. u,,(k) is a Bo()IE~all inpllt, e5e (k) is the output
of the event gem'rator, and YIJ(k) is the Bool('all output of till:' finit.e state machille~.
Thf' rules of the fillit,f' st,ate madlilH' can he:' l'l'pr<'sPlIted using thp logical relat.ionships
of SectioH 2.1. A t,\'pical finite' st.a.tc' Illachillr with rdgl's is ShOWll ill FigllI'l\ 2.2. fr tl\{,
s~:slc'lll is ill State' 1. it will transitiOlI to Sta.te' 2 as whl'n 15 1 = 1. It. will t.ransition
1H't.W('Pll States 1. 2, ami 3 aCT()nlill~ t.o the' valll('s of 15 1 61 , and (:1, hilt. State' 4 is
a tenninal state, 1>( c:ausl' onn' IIH' s:.,'st,PH! putprs St.atl' 4 th('I'p is no exit. transi tion
edge". This c!pfinit ion of a FS:\I is ('Olllprt'lH)llsiVC' ('!lough to ddiw' thC' l1lodels we' \';ill
llS(' later. For a mol'( gC'llC'fal, ft(\xihlr definitioll. Sl'f' [29].
2.2.2 Switched Affine Systems
GIll' form of all affiup c1,vllalllical systeTll is
x,.(k + 1)
Yr(k)
f(xJ.(k)) + ,(j(x1 (1.:))u,.(k)
I/(xr(k))
13
(2.10)
8 =1 2
State 2 State 3
o = 1
3
State 4
8 =0 3
o =0 2
State 1
8 = 1 1
o =Ov8 =1 4 4
Fignre 2.2: A fiIJit.~ :;tat.e machiue for a fourstate syStl~1l1 wit.h tat.(~ 4 a."l a t.el'l/liw:1.I stat(',
14
while a . wit h d affine y t In b r pr"",,\nb.·\d u in 11 £ 'Ll1 [ ]
Ai(k)Xr(k) + Bi(k)Ur(k) + B5,i(k)
Ci{k)Xr(k) + Di(k)Ur(k) + D;:"i(k) (2.1.1)
,,~here the set of linear matrice {Ai Bi , B ,i' i Di D5,d ha. the, ell" or i E
I = {1., 2, ... I}. \\ hilE (2.11) ha: often heen used t approximc t,(' the b havior f
the llonlim'ar SJ stem (2.10), based on some appropriatf' wit hing ruk wp will. h w
ill Chapt rs 3 and 5 how to H. e th . witched affine sy tern to mod I a .). tf'm " ho. E'
Iiupar d~'lul.mi(;s are affectcel by the state of th(' finit(' statE' machinf', and vice ver. a.
2.2.3 Event Generator
Thl' outvut 8(k) of an event. generator is a biliary sigllal that is g('nrratcd based Oil
sonw lille<ll' affine wllstraillt. for examplp, a threshold of statt' variahl('s x,.(k). Th~~
sigwtl 8(k) is all auxiliary variable used as an iUJ.mt to thl' fillitp stat.r. rnal'hillP. The
(JlItput. of the ('W~IJt g(IIl(ll'at.or is a t'lIrlction of t.he [(Ial signals in tlJ(' . witdl(ld rtffinr.
s:vstPlll 2.11. thus
(2.12)
TIH' Iwha"ior of all C'VNlt. g('Ilerat.or can Iw modeled IIsing til(' form df'scril)(\d ill
Sect.ion 3.2.1.
2.2.4 Mode Selector
Con\,prs<".v, t.he Ollt.put i(k) of the inod(' sdpctor is lIsC'c! to cont.rol th(' b('!Javior of' (.be
switched affine' syst<~1lI hasf'd OIl the lwha,\'iOl' of tIl<' finit(~ stat.(· lIIi:\('hill(', til(' values
of' xdk) and/or ub(k). or t!H' ('vent g(,I1<'rator. TllC'rc>f'ore, WC' hi'lVf'
(2.13)
L5
The mode selector awl it interaction ,;c. i h the
2.3 Chapter Summary
. ,ti lL .2. .
In thi' chapt r w 1m prest>nted au introduction to th l'l11('~ of Boo}pa.JI logi.c, Hnd
shown how arhitrar)' logica.l rules <.all 1)(' l'('prcscuted using fundameIltal relationships
as building block... W<" have provided Cl dt>finitioll of a hybrid s,vstern ill discrete h:vhrid
automaton form, and hown hm' its 'olllHOlwnt <,\1' oun rtf' I. iu . tllis form, WE'
ca.n capturE' a wide variet.v of Syst('Ill !Jl'haviol's. These ('on<:l'pts will allow liS to huild
detailed models in later chapters.
1G
Chapter 3
Hybrid Systems in Mixed Logical
Dynamical Form
We will hpgin \:vith a d(~s("J'j ption of how we ran east an optimization probkm over a
1l0lH'OIlWX s('arc.h spaC'P into a COIlvex spitce using auxiliary variable<. This provides
an rxample of {\ simpl!:' Syst!:'lll having rc~al state variaLle. alld hillary auxiliary variabl('
s, and the optimization prohlem provides an ('xi:uuplC' t IJr\J WP ('all solV<' ill ,hapt,('1'
l to <!('monst.ral.p stI'Ht,('gi('s for solvill!!; I\lIP problems. We' will t.IH'1I silO\\' JlllW DB
S.\·stPllIS call 1)(' n'p1'<'s('11 t.ed in mi:t:ed [ogiml dyuam'tcal (1\'1 LD) forlll. ast.i Ilg s.vs1.l'ms
ill !\lLD form [,('sults in a <:oll('is<' rrprc's('ntatiOIl t1lat VlI(' will ('xploit during scardl
algorithllls 1() find fr a....ihl(' solll tious awl sol V<' optirnizatioll programs. A bOllllcing
hall IH'<J\'i<lt's a silllpl(' (~xampl(l uf a dYIIHmic hybrid Syntelll t.hat we will d('snjh(' ill
l\ILD form and sol \'(' in Chaptc'r 4.
\\'e \\,ill lIS(' t he' program HYSDEL [30]. of['prcd as [1'('(' software hy pngill(I('rs at.
Il)(' Swiss Fpdl'rHI Illstitlltc of Tc'chllology. as all (lid in cast,illg h~rbrid s~rst(,IllS into
?\1LD form. Th(, lLD fonu cau j llcorporatc ('onti TlUOIlS dynamics all t,o11Jat.a . and
t h(' iut('l'Clet iOll h<'l\\'c'(,11 the two. Thc' rC'latioI1sltips d(~s('I'ihing MLD form allow us to
ulHkrstaud I h(' illpUt. alld output. of HYSDEL. \Vr will pro\'ide explicit calculat.ions
17
howiug OUI' boullcin . ball 'ample iu tILD £, I'm and h w h. H ,. L iupu ad
output. This will allow us t under. t.and how w call u HY 0 L to mod l the III r
rornplicat!~d RoboFJag ..irnulatioll robot l'oordillation control lJwbkm.
3.1 Convex Regions in State Space
For somE' real vector X r , jet. P 1)(' the dosrd, compact set of all possibl<' valucs X" can
takf:>. DrfiJlp SOIne partition of the, pace into 1 distinct convex rf'gions
(3,1)
as iu Figure 3.1 for I = 4. Let f(x,,) for X 1' E P be a piE'cewisl' liuear fuurtion so
that d.\'lJamics X,. = .AI X,. arc associated with region PI, x1• = 'lX,. an' associatf'd
'with region P'l.' et seq.
A COUV('X I'l'gion is one where ('wry point on a liIlt' C:Ofluectiug two poillts in the
r(;~gi()Il is also iu the regioll: so that for eV('r.v X 1',1 and X r ;2 in the' couvex region, the
poinl XL:I defiJl(~d by
(3,2)
is also ill the' ('onn~x regiull. For ('xampl,~, rWI)' [l'gioll ill FigllI'(' 3.1 is ('ouwx. '0.'(.
ma\' CtHltrast this with the l'<'gioIJ ddiIH'd h,v t hI' disj unctioll
(3.3)
for X r [ ]7' 'l []"I' (]7' = .1'1',1 ;1:1','2 E lR and a I = 1 0 ,(l,'l = 0 1 , bJ = 0.5 1)2 = O.s, as SllOWIl
in Fignre 3.2.
If I l)(' rq.?;ious are all rOIlVC'X, as shown in Figure 3.1, t.ht'n ('H('ll )'('gioll is dC'liIH'd
by H {"uuj 1II1('tiOll of tlie fOfIll
(3.4)
1
p
Figllr<' 3.1: A jJolyhedral part.ition of state spac(~ P an.:ordillg to linear tlu·esholds.
2
1.5
0.5
o
Xl
XI sO.5
X
2
$ 0.5
X2
0.5
1
1 0.5 o 0.5 1.5 2
Figur(~ 3.2: A disjUlldioll of t.hrcslll)lds ddillill~ a 1100JC()I1V(~X wgioll.
10
Each bOlllldary i E {I, 2 ... , I} i a h. (Jerplan z iat 'with a lin
(3.0 )
wjl h ~ollle eLi E }R1I aud bi E R or linear optimization algorithms s "lrching wi hin
a regioll Pi clrfinrd by (3.4), t.he general form of thE' program is
n~W1:IIL'tZ(:'
with Tes])cd to
.'j1J,ch that
X"
aTl X r  bl ::; 0
a7' 2 X r  V2 ::; 0
(3.6)
and we call seE' t.hat t.hrre is an implicit 'and' hrt\H'ell ("adl constraint pquat.ioll.
program l'C'pn!seIlu'd hy (:).G) is call 1)(' solwo by using oUP of Uli1I1Y algorithms that
('an haudle linear, quadrat.ic, or nonlinear objective fUllctiollS with liural' cOllstraints.
However, for systems that Illay indude wmCOIlW!X search spaces we cau solv' the
probl(l111 by drhuillg auxilim.v variables to cast til£' s. stem into a COUV(!X space iu which
we cau {'( H1V('lli('lltly coud uet au upt illl ihatioll IIsing a simpkx litH',u' prognmlluillg
algurithm.
If a tillite maxillluIIl or miuiulltJIl exists fur a liIl<!ar hUldiou f(x) wit II X r E 1R1I
withiu somr ('ouvrx regioJl P c!rtilwd by linear thresholds, th('n it Il111St. (~xist. at. some
H'rtpx of P [31]. Tvpical solutioll algorit.hms for linral' pr(Jgrammiug rpqlli['(! t.hr
s('arch rq~ioll to ))(' nJnvPX, howf'ver tht· shad('d regioll iu Figul'(' 3.:2 ddincd hy (3.3)
is dearly 1l0IH:on\'('x. By using an analysis ItlPthod similar to that disCllSSt J below,
WE' can cast a s~·st.(,H1 ddined by til(' combinat.ion of Ilf'uristic mit'S wit.h COIltiUU()IlS
('quat ions iuto a [unll ()(·fil1(·(1 by in('qmtlil i('s o!)('ratillg on I'I'a1 awl hiuar:,! variahl(·s.
TIl(' r<'lax('d roml ofr his r('pl'('sf~lIt.atioll product's a ('OllV('X statp spa('(' r<'gioB ill a
higlH'1' dinl<'usioll. and t.h(' l'<'slllting program IUlving hath n'nl al1<l hinar.v variabl('s
20
is called a mix dinteg ,1' pmgram ( iIIP). will r . urn 0 hi. x llpl. I ter in hi
chapter and in hapter 4 to illu 'tra e how a UP d fin .d in thi. rf'gl 11 c· n b .01 d.
3.2 Mixed Logical Dynamical Form
\ V<, call represent a hyhriu dynamical Systplll ill mi3:pd logical dynmm:eal (t.ILD)
form [32]. This form is couducive to solving larg(·'scak d~'lliuui('al s~'st('llIS, aull it
gives a s~rstel1lati(; way of repr 'elltin r til . rdati llships lWLWPl:'lI tIl(' ('OlllpOl}('nt.s of H
DRA clescrihed ill Figll[,(' 2.1. Th" form is dl'veloped using the logical and thn'shold
rPlatiollshi ps allloug tht· stat.e vectors X r and Xb, and t.he allxil iary \,(,('tors «5 and z
that W('1'(, df'vclopc'd in Chapter 2, III ord('r t.o transform a s.\'stelll in to 1LD fonu,
WP ul:'"d to d('hw' sewral rdationsllips in a cOllvClIiput form.
3.2.1 AnalogtoDigital Relationships
Th(' analogl'odigital relationships showII hC'I'(' can lw Ils('d to ddiJlP t 1[(' ('Wilt g('ucrator.
i,e.. tllPY dc\snilH' th(' way ill whic'h til<' sLat.(' of t.lll' swit,('!Jl'd <tffill(' S,VSI,l'11l ;tH'(,{'!.s
th(' hlli!(' stat(' nlil('hil1('. FllrUH'rmorc, W(' lila," Dud it 11,,('c'ssa.r~' to (')'(\at<' lI:n:l:ilim71
7)(l1"i.<LI)le.~ to n'pr0s('nt. t.h(' statl' of the s~rst."1I1. For ('xalllpl(' \\'(' lila.\' associa!.<' a Ji1.('ral
L with SOlll<' t lm'shold on a state' V(~ctor x. Sf) that
'( L = t1'l/(, H (I. X ::; I) (3.7)
",lwI'<' rr E JR1I is a v('('tor of ('odfici('lIts, x E ]R1I is II ,,('('tor of n~alvalll('d st.a1.('
variables. alld Ii E ]R is cl thresllOld .. If L bas all auxiliary variahiP r5. t.h(' ifandoldyil'
formulat.iol1 or (3.7) is l'quivalellt to till' iu('qualiti('s [30]
aTxb < 111(16)
(fT X  b > f + (11/  ()6
21
(3.8)
(.'1.9)
whflre 111 and m are realistic upper and lower bound of aTx  b' nd i. '\, :lIla.ll
tolerallC"fl, so that the r gioll II arX  bll < f i' negligible. (' can I'platp the. tatC:'
of x t.o the logical :t.ates by d .fining l additional auxiliar~' binary logical variabl('
(\ E {O, I}, i E {I 2:, .. , I}. A· we will how in Chapter ! b.v denllillg the appropriate'
auxiliary \"ariables W~ C,Ul cast the system into a form that call be ('ffectively !LaBelled
by MIP algorithllls.
3.2.2 DigitaltoAnalog Relationships
wllich t.he state of the fiUlt(' statE' machim. afh'cts the d~'muni('s of the switched affiuc
systelll. For example, certain r lation hips may be nefded in order to dpfine thf', tat.
if (S = 1 then (3.10)
the hybrid sYstem. ""Vfe' can shuw that in geBeral [30]
Digit.altoaualog relationships defiIlf' the mode, elector, i. P., thc~' dcscrilH:' the way in
of Romp auxiliary realvalued variable z E z based 011 t lteJ values ur 6 anel/or Xb for
is ('qllivalt'nt. t.o
(m2 1IJ1)() + z
(Tn,  M2)cl  Z
('11/ I  ]\[2)(1  r5) + z
(m"L  lIid(l  b)  z
.,. < (1"2 X,.  /)'2
< (J,T'2x,. +U'2
< "I (ll X r  b,
< a,7' X r + hi (3.11 )
where :~ is thl:' auxiliary n~cd variable, awl MI m,_ M'2' (\111/ 111'2 art' rmlistic llpp'r
ami lowf!r bounds for t 11(' pxpressions 'J" • 'j" (J I X"  iii and 0"2 x,.  U'2'
3.2.3 Logical Relationships
It may 1)(' that We' llPpd to define' logical rplatiollships IwtW('('Il j,illar.v nuial>lC's. \V(!
nlll llSP thE' rdat.iouslJips (l('fin('d ill Cltapt.f'l" 2 as building blocks to cI('fiuC' all.\' logical
22
infC'rerKr machinf'.
3.2.4 Continuous Dynamics
\,Vt' ClUI mod(') th(' hrhaviOl' or thC' A hy df>fillillg the dynamics or a . ysterJ) in t.erms
of the i'luxiliar variables ill addition t.o th tat.e and input. .ctOl'S, 'md lLsing t.h{'
~ I to (,oIltroJ t IH' allxiliary ,'aria blE's, The eli 'cretized st.. .spac<' repres('lltat.ioll of <:l
s:vst(~111 ,,\Till havc t.he gell('raJ form
x(k + 1)
y(k)
A.x(k) + B[u(k) + B16(k) + B3z(k)
C'x(k) + D,u(k) + D;l6(k) + D:1z(k) (3.12)
whC:'re u(k) [u r Ub] is a vector of iuputs allu y(k) [y,. Yu] is a v ctor of
outpnts.
3.2.5 Mixed Logical Dynamical Form
III addition t.o the' rdatiolls ill (2.9), (3. ). (3.11), anu (3.12). W(' can rtlso drfinr lillC·rtr
a.Ig('l>raic l'rlatiollship. alllollp; <'1 H1PJlts of X,., linC'ar const.raint.s Oil Xl" and logi('al
('onstl'aillts 011 XIi a.nd c5 ill a fairl\' st.raightfonvard ma.llIlC'r. \V(, ('atl ('())lIhilll' th('s('
rC']atiollships to silo",' tht· s.vst('l1l ill ?...lLD f0['111
1.;(/;+1)
y(k)
E;lc5(k) + E:3z(l.. )
Ax(k) + B1u(k) + D;lc5(k) + D:1z(k) + Dr.
C'x(k) + D1u(k) + D;lc5(k) + D:lz(k) + Dr.
< E[u(k) + E.1x(k) + E"
The ('(·latioHships giH'lI ill SI'('tiollS 3,2.1 to ~L2.3 (','tIl h(' ('olllhilll'c! with all~' ot.h('r
approprlat.(· (,o!lstraiut.s t.o fOI'lIJ t11<' ill('qllaJit~· ill (3.13), wit II I'Odfi('i('lIt.s giVl'lJ by
th(' Ei llIrltI'i('(~s. Thl' offs<'t,s D" alld D·, do ltot iucn'asr the ahilit..v to l'xpn'ss S.VS!('1tI
l)(~h<l\'iur. IJIlt a/"(' illc')ucl('d h~' [28] for (,OllV0.Hic'!I('(', TJ)(, llot.a.tiOll Ils('<! h~' rbI' HYSDEL
output allows us to U10I'(' Ill'('('isd., pxprC'ss the structUI'f' of tbe 1l1()()l'I[28]:
23
[
y,.(k) ]
Y/i(k)
( .1)
""here the yahlPs of thp matrices cH(' oefiued by t II(' previous E'quations. s we ,hall
see ill lat('r chapt('l's, typical s.vst.ems can have' t.hollsallds of ('onstraints and variables,
so that pllttilll!; the S~·St.f'1l1 in this form and fimJillg til(' solution call 11(> a daulltillg
task. The program HYSDEL nlll IH' lls('d to transforllJ ,I s.vstel1l frOJll a fairly illt.uitivr
spt of rplatiollships (() t be ,1LD fonu of (3.11). B.y tasting om syst(~IJ1 into this sim pl('
form: we' makr it lILOH' allH'na.blP to solving state'spacp Sf'ardlillg awl optilllizatioll
algori thlllS, as W(' slIa II s('(' in ChaptC'rs 4 alld ;j. First, W(' will pn's('/Il "}Laud ('al
(·.IlatiollS'· for two Si/llpl(' ('xaIllpl(~s to ilb/slrat.(' t.lll' pw('('ss. In Chapter ..\.. WI' will
(~xa/Ili/)(' t})(' HYSDEL i"pnt. and (J1ltlJllt fill'S ill iliOn' dl't.ail.
3.3 Example 1: An Optimization Over a NonConvex
Space
If WI' wish to maximizt' an arhitrary flllJ('tioll of two Yariahlf's :1:1 and ;r:'L ()V('l' a nOIl
('OIlV('X span' such as that i.n Figu[(' 3.2 WP Iw('d I' 0 r<'st.at(· (3.3) usi "12: tlL(, alla]og
I'odigit.al a"d logical rl'1at.ioTlships d(~s(Till('d ahoVl'. 'otC' that this is noi iL dnlHIIIi('
syst('lU, so \H' ('cUllJot call it a mixf'cl logical dJj1l.rLrnir:al ('on)]. How(·\·('!'. it is a good
illllstrat.iorr of t.1H' logical reprf'sentation, and p;in's a simplE" ('xampll' to illllstratf' t.lH'
t ('dmique' of soh'iug mixed int!'gp)' jJrop;ralJJs as d('sni l)1'd iII the' llPxt chaptN.
The di jUllctioll d ribing h ff'a ible (,Cl. J
. ( ,15)
for X r E ]R2. D('finl' the auxilia.ry ,'arialJ1<'s
6, 1 H (/ I X r ~ II I
151 1 H (L1XI' ~ b2
():\ 1 H (fIX,. ~ bl V (/2X" ~ b2
(3.1G)
with (/ I = [1 0]. ([2
equatioIls
2 . This can \>(' rq)['(~seuted by tIll'
niX,.  bl < H:r1 (161)
U1Xr  U2 < l\lx2 (1  61 )
(JlxrIJI > f + (m.d  f)6,
(L2Xl'  U2 > f + (mr'l.  f)62
A'I + 62  6:\ > a
()I  6:\ < ()
62  rl:l < 0
():I > 1
( .17)
(3.1 )
( ,1 a)
\Yp ha\'(' HO dynallli('s for tltis ::;ystem. so we call l'l'}>n's('nt. tiL(' S,\'stc'ltl ill t.II(' fOI'l1l or
(3.14) as
Jel () () 1 0 IJI + MJ"
0 .'\1.,.2 () (J 1 112 + [.1'2
'" r l  f II () 1 0 hi  (
() Tn.,.2  t O
t5~
(J 1 h2  ('
(3,20)
1 1 1 0 ()
X+ 0
1 0 1 0 () 0
() 1 1 () (J 0
0 (J 1 (J (J 1
\\'Lll'l'{' 11('1'(' X,. = [.r I .l':,J'. and 15 = [6 1 (51 6ll', rr tllf' binar,\' variables are' 'rl'iax('<!'
to cOlltimI<.Jus variabl(~s 15, E [0,1 fJ. \VC' call dc'fiIl(' a !lew colltill1I<.JUS state' variab)('
25
v =; [X,. Or] E jR+2 X [0,1]3. TIlell till' inequality in ( .20) c' n b(' writt. n .
50[[1(" Av ~ h. Thr lI('xt cha.pter will show how an optimizatioll al . riUun can be
implemented in this space'.
3.4 Example 2: Modeling a Bouncing Ball
III order to dplIIollst.rat.!' one approach to modeling a hybrid s.vst(:'tn W(' 'will pl'(,s<'llt
th<' hOUllCillg ball exa.mplp. This is a C0lJ11110n t.ut.orial ex;u11ple of a hybrid systelll
[29], [30].
A hall dropped above a fiat rigid smfacC' will fall onto the slIrrae(' a.nd UOUllC('
Oil impact a.ccording to sOllie coefficieut. of dissipation. This is a h urid dynami 'al
systelll with the lOntinuous state space elivided int.o two distin t n'gious. When the
I)aU is not. in contact with the surface, its motion is modeled hy
h v
(3.21 )
wlll'[,(' Ii is til<' hpi~ht of the hall ·tbo\'(' t.lI(' surfaCl', 'CI is tIl(! Vl'locil'.,\' of t.lw I>all. III is
t hI' mass of th!' hall, awl .II is the ;u'('deratiou COllstallt for gravi t.v. TI Ie approxilllfl t(,
solutioll t'O 1h(' clis('I'c,t<,tiul(' equat.ion at time k: is
h(k + 1)
v(k + 1)
h(k) +1:v(I.')  T,';.II
v(k)  T.~(J (3.22)
""'hc'l'l' T" is a timc' interval for ('omplllatiou. UpOlI impact, tlH' ball ['rve['sns dir('ctioll,
and <) dim'rent d,VllflUlic eqllatioI! des<Tibes t.hp s.VSt,('1l1 c1l1l'ing the' illlJ>(l(,l. tim(' int.l'rval:
h(k + 1)
v(k + 1)
Ii(k)  T"v(k)
(1  o)I'(k)  T,.q (3.23)
[h(k) 'dk)]"', t.hpI! we han' olll~· OIH' linpar thrps!Jold couditioll. I,(,(Jlliriug
2G
()11~ l>illar~' variablf' 6(k), such that
11np(Lrt . t'l"LL ++ 6 = 1
6 = 1 ++ [1 O]x1·(k) ~ 0
\~ hich is <'qui valeut to
[1 O]x(k) < h1/UtAl  c5(k))
[1 O]x(k) > f + (I1 171 i1l  f)(5(k)
111 this ca~(\ W(' call defillc z(k) = [h(k + 1) v(k + lW', then
(3.24)
( .2'")
(II/',!  M,)6(k) + z(k) < [1 T'] Xl'(k) + [ +1'';.'1 ]
() 1 +T,y
(11/1  III2 )fJ(k)  z(k) < [~/ __;'] x 1·(k) + [ it:; ]
( III I  AI'1 )( 1  c5 (k )) + z (k ) < [~ _(~ ~' n) ] x,. (k)  [ :~ ]
(1112  i\/, )6(k)  z(k) < [~)l (t:;t)] x,.(k) + [ ~~ ] (3.27)
if c5 = 1
is equivalent to
flle'l/,
el:w
z(k) = [ ~
z(k) = [ ~ (3.26)
wll('I'('
1//1
(3.28)
, otin' that (3.27) prl'risC'ly derrllf>s (as all (\CJuillit~,) t.hl" "alup of z(k) "vhf'lIeV(~r
6(k) is fixed. If ~(k) = 0, t.hen tht' first alld s('colld ill('qualiti('s form au equalit..,·,
and t.he third and fourth equalities arr trivially true. C()lIvers(')~, if (5(k) = 1, th 'U
t.he third alld fourth iIlf'quaLit.ics form a.u ("quality, ilnd the first and s('("oud equalities
~trr triviall" tnH'. Vve know this must ]H' so, because the hOllllcing ball has oILl,v OIl('
feasible solution for a giveu S(,t of paranH'tcrs and initial ('omlit.ions,
For our !>ollnciug haJl example, the disrrrt.rtim(\ rl~·na.ll1ics alld t.hr rclat.iollshi ps
iII (3.25) and (3.27) d('fin(' tile !v[LD form
2
h71IrLT
(h7/1i1l  ()
h111i1l + T..7'1IIill  T}',fJ + T.,'Umill  "'7/1fu
7''l/Lill  T.~!J + (1  O)'Cm i1l
1I111i1l  T., 1'"IU.,.  huwr  T"UULlLJ' + 1'/',lj
(1 n)V7I/Ct.T  "111el:!' + T,y
h1l1i11 + T,ullUL.,. + htrla ,,. + T"I'/ltllJ'  1',;y
(1  n )'U/TLfI.,. + l'rnar  T.,g
h."'in  T.,t'fIIi" +7';y  T.,'llll.ill + lI met:]'
(111111 + 7~,q  (1  U)'/m;"
r5(k) +
o a
o 0
I ()
o 1.
1 0
o 1
I 0
o 1
1 0
[) I
z(k) :S
n 0
1. ()
1 7,
0 1
1 T.. x,,(k) + () 1
1 T,
0 (ln)
1 T,
0 (1  n)
hl/w..J:
(
T,'2Y
1',y
T,.'2IJ
T.,q
hmin + T..i:i'IJ7I/.(I.:!' + hllLa./, + TS'(lIULJ'  T/'g
(J  n) 1J7l!<13' + 7 111u:/,  T.,g
!l.1II;in  T..'l'miu + T}g  T.,'17min + hl/luJ'
Vmin + T.,g  (1  0')'( min
(3.29)
The solutioll 10 this system is unique, that. is, giyt~Il vah.!ps of x,.(k) and 6(k)
precisely ddiuf' x,.(k + 1) and 6(h: + 1), \\p will show ill thr lIext. ('hapt,pr !low S('iUTIJ
al/';ori t hillS fiud til<' solll t.iou t this Systl'Ul of (''Illations,
3.5 Chapter Summary
In this ciJap/'N \\"(' han' sho""u how to c'ast hybrid syst<'InS into llIixc'd logical dynamical
form, aud UPl1Iolistratc·d hm·,' this ['OrIll is a HC'xihJc> and powerfnl \\ay to r<'prc's(,llt thl'
1wlJa"ior of a ud i IItc'faction I)('tW('('H huit c'sta tc' lIladli Ill'S cUld switdl<·d aHiw' syst(IIllS,
:\lallY n~al S~'st(,lIIS nUl 1)(' (Ilodded ill t,his form t,o capture' cktflikd 1)('lJavio)'. illdlldillg
multiagt·.!)! rulHlti(" s~'stl'lllS with cOlllpl('x goals. . CkUIOlLstratioll of SlId. a syst('IIJ
is thl' sulJ.i('('l of Chapter G, Bdor(' WC' dn that. \\,(,'11 sho\\' IIO"" to solvc' IlliXl'diutc'gr'r
programs (lsillg hnull"halldbolll1cl lIH't ]J(Jds ill Chapter ..1.
29
Chapter 4
Solving Mixed Integer Programs
rIavillg d(.'hllPd a hyhrid SystPIll ilud shown how it. call IH' cast. illto '1LD fOrlH, we' !LOW
discuss lim\' the a mixedintt'gel' prognuH call 1)(' solved Ilsing till' bmnchundfwund
t.t dmiquf'. \I\'r in trod uct' t IH' Silllpkx algori r.hm for solving liw'ar programs. and shmv
how linear jH'O/SJ"{IUlS form sll1>prohlc"llIs for t.ll(~ ImulchaudIJolllld algorithlll. \V(,
show tlH'sr t<'chuiquf's primaril,\' by dpllJOIlst.ratioll, for which W(' \vill Rolv(' Example
frolll Chaptc']' 3 wit IJ (l lilH'ar ohjc'Ctiv(' t'ullctioll ov(']' a I1OJl('OllV x s('ardl spat'("
TIl(' blllllicing Iml] (>roh)(,Itl. Exampll' 2 from tlu"ltapl.('r 3. willlH' mwd t,o dl'IIIO/lstrate
till' plO("('dun' for soh'ing; r.lIP's IlSillg the C()I1I1lH'ITial sol\'('1' '1 LEX availal>]l'
froIll ILOG [33]. \\"(' \\'ill traIlsfofmillg it. into l\1LD fOl'lII lISillg HYSD  L, t.h('lI uSI'
a program written for this tlH'sis IIsing \IATL B soft,wan' t.o writl' t.I11' MIP inpllt.
fills, \V(' \rill solv<' tlw prohl(>nt llsing a CIrlllglJag(' progralll which calls t1H' CPL X
lihraric's, \y(' \\'ill S(>(' that we call so)\'(' t]1<' !lOllTlcillP; hall probll'JIl, wlticlJ is both
rrclll'siw (ill disrn'te tiuw) auc! couditioJlal witho1lt a sinp;ll' "for"' I<HlP or ·'ift.IH'Il"
st atc'lJ)C'ltt b,'" usi up; t]1<' i\ ILD forUt.
30
4.1 The Simplex Algorithm and Linear Program
. IDlng
A silIlph~ lillear program ov('r a cunvex r('gioll P with rf>sIH'ct to X r E ]R7I is \\'ritt('1l
111m?',Wl ;;('
with 'I'C'8fH'd to
.<;/tch tI/(d
(4.1 )
wlJ('n' f (xT·) is a Ii Ilear olJ.ipcti VE' functioll, A Rwl 13 C(1 Il ddill(' liIl<'ar thrpsholds and
..y C jRft is t1l(' sd of all values that. X.,. is allow('d to take. Thl' Illatric('s.4 E ]Rlx,!
and B E ]R1I defiul' a sC't of J /)olllldari('s of a ('OUV('X rC'gion. ",h('n' (~al'lt boundary
i E I {I, 2, .. " I} of the conn'x ['('giOll is a h.vp('rplam· associat.('d with a liuC'ar
tbrpshokl a/x,. S Vi, with SOUlP. {/.'i E ]R11 and h; E lR. Thc' eqllatiOlls ; of tl](' lillPiil'
thresholds form rows of thl::' matrices A and B.
\ire will dC'1lI0llstratl' thc' solution of this prohlelll h~r t hl' simple:r method. Th('
silllplex Illf't.hod wa.." c!C'wlop('d ill the 1940's and llladC' famous I).\' Danzig, wllo ntlkd
I,he' probl('lll it solves a linem' fJl'oym,m[31]. This U'rtI1 pn'dates tIl<' widespn'ad IlS('
of comput('r programs. anel it 11 II for t.ll II atcly has a similar uaulC'. WI' should silliply
lllH1l'rst.cwd the diff('n'I1(,(' I)<'tW(~('11 a lilwa r program and th(' COIl1 I'll (.c·r progralll that
mig/II' h(..I used to solve it., AIl (lxamplp of t.ll(' silllplc'x aJgorit.hlll appli('d t.o Ollr
optimization ov<'r a 1I01H'OllW'X span' will J)(, SlIOWll fw!o\\'. Hc'J'(', W(' will cksnilH'
Itow t h(' prohlclIl is cUllstrnct,('cl.
The' simplex algorithm is a syst('lIlatic BH't hod of sol\'illg hllc'at' pmp;ralllS ill tit<,
fortll (4.1) with ('quality or ill<'qual its (,()lIst.raillt.s. For ('\,pry i Il('q llal it\' ('OilS! ra iIll, wC'
will illtl'O<l11(,(, a slack vur'icLblc t.o C<H1\'('r1 it to all ('lfllalit,\' (·OIlstraint. Tc'xtlHJoks Oil
tIl<' suIJjc'rt t.~rpi('a.Jh· i.ntrodlln' awl d('lllulIst.rat(' tllC' lIJ('thod 1>~' \Va." of ('xillllplc'. and
"'(' "'ill do t 11<' saJllt', Th(, simpkx alp;ori t 11m ih ;\ nI<'l bod for ('(J!l\'('X opt illl iza t.ioll.
however, we are interested in op imizatioll r 11 lICOll x pa.. h r for
will demoll tm e his method a a subpr bl III of larger UP soluti lid' rib l ill
the next Ectioll.
4.2 BranchandBound Methods for Exanlple 1
AlIlOIlg COlIll)('t.ing lll<,thods for sol"iIlg !\np s most. CUll H1H·n.:i al software' packag('s
(illtludiJLg CPLEX) use SOIll(' type of unJ.nch.n.nduound algorit.hlll[27]. Th(' basic
iclPa is quite simpJr: first. t.lJl' "LP RdcL'{(itioJL" is form1lIatC'd, w1le1,(, all hiuarvntlu('d
\'al'iablcs an' allowed to take <tlI~' rcal vallI<' ill [0, 1]. Thc'I1 til(' optill1izatiou is
solved as if it were a rcal\"alllpc! LP problem with thc' appropriate b01lIlds on tlle
st.ate variables, indudiug thl' relaxed binary stat varialJks. If the n~sllit lraplJl'Il 'd to
show all hillaryvallled \'ariahles havilll!; integer valm's (pit.hl'r () or 1) at the optimum.
thr11 this is the glo] Jell opt ilIlUIll for tiH' s~rstelll and the algorithul 1,enn inatC's. Sill(,('
this is Ilulikdy ill ;\ liontrivial l~xalllple, tlw algorithm choos('s ow' of t.hr fractional
variables and neat.(·s two sllbpro]l]P!llS, OIl(' assllllling it is has il valllC' of 0, 1.11('
ot h('r assnnll'S it ha,s a ,,'altH' of 1. EaclJ S\l b probll'lll is sol \'('d \ alld flllthc'l' bralid ling
I)(TllrS as tbl' opLimizulion bn; is huilt. As morr variahll's ,Ul' fix(·d t.o int('gl'r valtll'S
awI furt hl'r SI1 h ]>l'Ohl('llIs an' soh'HI. t1)(' braudl('s of' 11](' opl.i III izat.io!l t,l'(,(, grow.
Feasihk SOIIlt.ioIlS an' t.hos(' tliaJ satisfy t/w cOllstraint.s with SOIll(' l'<'lax('c! valll<' of
t.he "ariables. Owr a f('itsi hk su[n t.iOIl wi t.h all hi liar." values is obI.aiw'd, there is 110
IJ(,(~<! to (lxplore a. givC'1l braudl ,lilY furt.h(·r, aud the branc·1J is said to 1)(' faf.lwTn,:ri. TII('
first .feasihl(' solution wit.h iutl'gel values oht.ained h(~conws tIl(' ·i.nr:mnlwul: solution.
which is H'tul'llrd as tltl ' OptilllllJrl if lat,l'r calculat.ion fail to SIIOW a fl'<lsihlr int.c'/.!,er
solution wilh <} Iwtt.rr obj('l'liVl' fllUct.i()ll. Any slI('rrssiVl' feasibll' il1tl')1;('r solutioll
with a h(,tt('l' ohjl'nin' fUlil'tioll \'aI1l(' will <!P!>OS(' thl' old iIlCtllIlI)(IIl1. solnti(!lJ. All
optilllizatioll tn\(' wilb llotat.ioI!s is incllldc'd with t.}l(' ('x<\lllpl{' solwd ]w]o\\,.
32
Th(' dC'fault setting of t!Je PLE T r..fIP I .r let. the pr gTam eh 0
of solution for the LP 'uhproolplIIs I>a.<;ed on its own crit('ria' l1s11all~' tIll' dual simple'x
llIethod is dlOSPli [33]. Our tri"ia1 (·:;:;:1,]1I.plp can he solved b~ ins\)('ctioll. hllt, in a I· rgel'
systCI1l efficient mpt!Jo<!s of findiug sollltiolIs will be V('r~' importaut.
In our f'xample, the' solution of each subproblelll is (,;lsil.v ohtailled IlS111g the'
Optimization Toolhox of IATLA13, or l1~illg a uumber of LP t.('dlllirl'l(·S available· ill
introrludorv textbooks OIL the sl1hject, \,Ve will show t.!J(' progn'ssion of I hI' soll1 tion
llsing thf' hranchandhOllnd t".<'dlJl iqllr. lea"ing 011 t the soll1 (iOIl of padl Sll h pl'OlIl('1I1
after t.ll(' first. This will hel p ('xpl il iII and ulOti va te om descri pI iOIl of t hp olllpll t of
thf' CPLEX trials showlI ill lat,l'r chapters.
4.2.1 Optimization of a Linear Function III a NonConvex
Space
For t.he problelll of Example 1 int.racIlI("I'd Chapt.er 3, define' v
[f we rplax rll, (2, and 15:1 to han' )"('HI vallll's ill [0. 1]. w(' obtain tlll'slItlprogl'alll
al ( II(' first node of tlt(, lllixl'dilltf'g('l' pro/1;ralll. To ('lllpltasiz(' tlt(, ['ad t11al all IIII'
\'C:tl'iahlc's at nocle. 1 (H(' t.n'al(·d t.1ll' sauH', l'x('ppl. tha t lowpr all( I II PI )('1' hOlluds or ()
and 1 an' applied to those' relax('d hillar~' variabl(~s, Vi(' will rpf<.l' t.o t.!I<' (·O!lIIHlI)(·tll.S or
v as v = [VI 1''1. .//:) VI 'i!;Jr. \""(' wiU c!('s<Tih(' how this illitial t'f'laxatioll at. !Iocle 1 nm
he' solv('d hy halle!. This will allow liS to illtrodu('(\ I('I"llIS awl ,·itriahl(·s, ;wcI Iwlp to
f'xplaiu SOU\(' of'rJH' steps the' CPLEX prugralll gO('S tllnmglJ wlwlJ it is solvillg a MIP.
This is important h('callse CPLEX allows a /1;T'<'at. dC'a] 0[' f!exilJilit,r ill Ilo\\' W(' iIlSt.l'IIC(
it 1.0 soln' ('ach problelIJ. I'll<' stawlanl verbosp CPLEX 01ltPllt bas 1I1cUl.," 1('1'1 liS and
n\lurs. awl t }1l'S(' ('Xilmp!<'s will Iidp liS h('gill t.o lllld('rst.aud t 1)(' hasil" asp('l't S or t he!ll.
4.2.2 The Simplex Method Applied to a SubProblenl
\\'r' call represent t!lP inpqllaliti(~s of equation (3.20) at node 1 by (,Ollll>illill~ t I\l' state
wctors, insertillg thp vahlf'.
II!:!, I 1;0'2 = 10
?/I,x l 'IIL J'2 = 0
bl 3
b2 2
aud rr'arrallgillg:
1 0 10 0 (J 13
0 1 0 10 (J 12
1 0 0.1 0 0 3.1
0 1 0 0.1 0 2.1
0 0 1 1 1
v<
0
0 0 1 0 1 0
0 () 0 1 1 0
0 0 0 0 1 1
(t. 2)
(1. 3)
when) tlJ(' last. equat.i()ll ewmres that 'U,. = 1. This is a s~'stpm wit.h ;j llukuoWIlS awl
8 illpqualitic's, ami TfJ71k(. ) = G. To fa(;ilit,atc~ t.IIC' sollltiCH! wc mllst int.mdll('l' slad.:
'IIflTiables whirh rq)l'('sput thc' diH't H'llC(' Iwt\wc'll til(' !l·nhalld sid(' auel righlhalld
sidr' of cal'll f'quaticlIl. For exampl(" tlH' first ('Cjllatioll IW('UllH'S
(1.4)
III adclitiull. we will J'Cprc'Sl'llt t h0. ohjc~d.i\'l' fUlIction vallll' as a variahle· 110 which C'HlI
1)(' illC'orporatc·'u ill tl1(' syst(~1ll of' eqllations. \VI' will switrh to ('qllatioll forlJl ill onkr
to hdp llS kec'!> tra,ck of c'arh vari ablc' for dis(,llssion.
First. sincc' \\'(' kl1o\\' that the pighth ('qllalillll JIIC'rf'!:v sC!ts 'II;, = L w(' will c·lilllina(l'
/1,, fJ'Olll t11c s.vstC'fll. This is ('1111('([ tltp pn:solvc; stage ill opt.illlizat iOIl parlauc'(·. aIHI
CPLE.\. always c1m's this WhC'Il('Wl possilJlf' to H'c!lH'(' tIl(' sizc' of t.llP Illatric'c's ill ('act.
it ('aUIlOt. 1)(' t.otall~· disablc'cJ 1).\' thc' IlsC'r. Thc) n'sllit is
34
Va +1211, +2u:l
81 13 VI Hh:1
'';2 1:2 V2 10111
,';:{ 3.1 1'1 0. 103 (4.5)
.';4 2.1 (1"2 0.10"
8 u 1 1':1
87 1 'I'~
The variabks OIl t h<, LHS of t.1l(' ('quatiOll ,U'(' caJl('d !Jo.sir v(£1'iablf'.,'; whill' thos<, 011
thE' ri,ght are ealled rW11,Vasit: var·iablc~8. oilbasic \"ariabks must. 1)(' at l'itlH'l' Ht
their lower or IIp])('r hallUc! hel'e' th('\" an' all at t.heir lo\over bound of Z('l'O. Fwm
this poiut, we consid('r tll(, slack \'ariables and objC'ct.ivc function Y(\l'iahlc' to 1)(' .iust
like' thl' ot!lcrs, with till' ("ollst.raillt tllat 8; > () \;f i = {J." .. 7}. 'vVr> flOW haw a
svstpm with 7 c'Cjuatiolls and 11 Ilukllowns IllCl,kin,g it ulIdC'rdrt,erInill('d and elllowing
liS ·t c1cgn~l's of fnwlolllill optiruiziug ou!' S~Tst('UI. \V<, ("onsider it trial SOIU!ioll as
OIH' wlwrC' all til(' !loubasic: variablc's arr zero, and thc' basic" OHf'S an' 1I0UZl'l'O. This
solu hall would 1)(' h~asi hi" c'xc'('pt that .'i:j, aud .'i"1 an' llrgati \"f'. \V(' haw to introd uC'p
a S(lt of aTtz./ic'iul '/J(LTwbles that repn's('nt t.llt' diffrrcu('c·' brt.w(,pll t!H' hasiC' variabk
va)lH' (Iud i'\ fC'ilsih)(' nthll'[3':)]. HC\H'. IN t.hc\ hasiC' \"ariahl('s 1)(' <1\ = ()  .'i" a\l(1 \\'C'
will millilllizc tlI(' SHUI ofinfl'asihiliti('s. = <P;\+<PI' Silll'c" fot' C'Xitillplc', .';:\ = :J.t.
wp haw <P:I = 3.1 alld 8:\ IWCOIIlI'S ,I nOIlhasiC' vilriahle. Sifllilarh', W(' fllnkc' 8 I a
u(j[Jlm.'iil· variahl<'. so I hat S = '';:1  8.1. alld
'Po +12uJ +2'/}"2
8) 13 I') l(h 1
82 I:.? 1''2 l(l'!,'1
<JLI :1.1 III (J.l'P1 +8:\ q) J (t.G) 2.1 I''L 0.1/'1 +·'il
8(. 1 11;\
.'iT 1 1'·1
':),:2 'PI 11"2 0.1 P:j (J.11.1 + ,';:\ +'';1
Iu p}w,'ie one of t ht' soilltion \\"c. willlUillilllizl' so 1hal thc' artifiC'ial ntri;lhl(·s arc' all
zero, at which puint w(' will !la\'(o a fl'(lsil>l(' solutio!l. \\'p will illC'!'c'i\.Se '/'1 as mucb as
possible, so t hilt cD:l fl, alld
30
Va 36 I2<P3 + 2V2 +1.2V4 + 1283
81 a +<1>3 +9.9v4  3
82 12 V2  lOV4
VI 3 <1>3 +0. IV4 + 3
<P4 2.1 V2 O. Iv4 +S4 (4.7)
V3 1 'V4
86 a +V4
87 1 V4
S 2.1 + <I>3 V2 O.Iv,! +84.
\Ve will increase V2 until <1>4 becomes zero, so that
Va 40.2 I2<I>3 2<1>4 +V4 +12s3 + 284
81 a +<fl3 +9.9v4 83
82 9.9 + <I>4 9.9v4 S4
VI 3 <1>3 +0. IV4 +83
V2 2.1 <I>4 0. IV4 +84 (4.8)
V3 1 v4
86 +V4
87 1 V4
S +<1>3 +<1>4
ow S is zero and (4.8) represents a feasible solution with objective function value
40.2, Xl  3, X2  2.1, c5 1  1, and 62  O. This is th end of phas on I
which shows that a feasible solution exist . It is not optimal, since our goal is to
maximize Va and we could increase it by incr a ing any value in the xpr ssioll for
Va that has a positive coefficient. This begin phase two of th proce sfinding th
optimal solution. We could, for example, increase s::\ until Sj becomes zero, so that
Va 40.2 2<1>4 + 119.8v4  1281 + 284
83 0 + <I>3 +9.9V4 s]
82 9.9 +<1>4 9.9v4 84
Vj 3 + lOv4 S,
V2 2.1 <I>4 O. IV4 +84 (4.9)
V3 1 V4
86 0 +V4
87 1 V4
S +<P3 +<1'4
VVe can increase V4 until VI  10. Recall that a nonbasic variable is either at its
lower or upper limit. In thi case, our logical formulation has an upper J)()und of 13
36
but the actual bounds on VI
at a value of 10 to yield
Xl are 0 < VI < 10. SO VI nter th nonb t
Vo 124.06 2<1>4 0.02s] +2 4
s3 6.87 +<1>3 O.Ols]
S2 2.97 +<1>4 0.99s1 S4
V4 0.7 +0.1.'1]
V2 2.03 <1>4 O.Ols) +S4 (4.10)
V3 0.3 0.1 )
S6 0.7 +O.ls]
S7 0.3 O.ls]
S +<1>.3 +<1>4
Finally, we can increase S4 until S2 = 0, so that
Vo 130 28)  282
S3 6.87 + <I>a O.Ols]
S4 2.97 +<1>4 0.99s)  '2
V4 0.7 +O.ls)
V2 5 s) 82 (4.11 )
Va 0.3 O.ls]
86 0.7 +0.ls1
87 0.3 0.ls1
S +<1>3 +<1>4
This is feasible becau e all of the variables take values within their allowable rang.
Notice that both of the nonbasic variabl.es in Vo have negative coefficients. Increasing
either of the nonbasic variable would lead to a decrease in the value of the objective
function. Therefore, v = [10 2 0.3 0.7 l.Or with an objectiv value of f(v) = 130
is a feasible optimal solution of the relaxed problem.
4.2.3 Exploring the Optimization Tree
In Section 4.2.2 we have found a feasible optimal. olution to th relaxed problem,
but not a feasible integer solution. It is shown as node 1 on the optimization tree in
Figure 4.1. ote that the value of the objective function for this trial becomes all
upper bound for the feasible integervalued solution. This is because when w restrict
the binaryvalued variables to carry only integer values, the solution objective function
37
Node 1
z = [10 5 0.3 0.71]'
f(z) = 130
() =a 1
Node 2
z=[102011]'
f(z) = 124
Feasible integer solution
becomes the incumbent
, Node 3
z=[310 10.121]'
f(z) =56
o =a 2 () = 2
Node 4
z = [3 10 1 0 1]'
f(z) = 56
Feasible integer solution
with suboptimal f(z)
Node 5
z =[3 2 1 1 1]'
f(z) =40
Feasible integer solution
with suboptimal f(z)
We will choose to branch on 61, creating two Sll hproblems. Th first has oJ = a
Figure 4.1: Optimization tree showing the decomposition for example 1 using the branchand
bound technique.
subprogram using the function "linprog" from the MATLAB Optimization Toolb x,
\vhich solves the prohlem of equation 4.1.
[10 2 0.0 1.0 l.ojT with j(v) 124. This is a
must b Ie s than or equal to that of the relaxed problem. ow that we hav se n on
method for solving each relaxed subprograms, we will solve each of the r maining LP
and yields node 2 with v
feasible integer solution that becomes the incumbent. solution. 0 further testing of
this branch is requiredwe need not test 62 in this branch because a feasible integer
solution was already found. This branch has been fathomed. At this point the gap is
the difference between the upper bound found in nocle 1 and this result, so the gap
10
,'.
........
",
, ' .' ,'....
.." ' '.
, ,
~ ", fIlode 2
. ~ ..
".:
....
.': .
.,:.
o 0
·... , .'
140
':
~ 120
)( ", :
;;: "
,.'
100 ....
"
80
60
" ,.'
40
20
0
10
Figure 4.2: Objective function value and bounds for example 1. Feasible solutions at various
nodes are marked.
is 130  124 = 6. Th ratio of the gap to th) absolut value of the solution i often
used as a. convergcnc crit rioll) and it is part of tlJ' report d output of PLEX. The
solution for node 2 is shown in , ib'1U 4.2.
The s cond subprogTam has 6) and yields v 13 10 L.O 0.12 I.Or
with f(v) = 56. This is not a feasible integer solution, sine 62 does not llave all
integer vdlue in the solution to the relaxed problem. We bra.nch on <52 while ke ping
1 ~ required for this branch of the tree. Th olu ion with <52 oyield
v [3 10 1.0 0.0 LOrI' with J(v) = 56, and the solution wi h 02 I yields
v (32 1.0 1.0 l.O]T with f(v) = 40 which is a feasible integer solution to the
optimization, but has an objective function value less than our incurnb TIt solution,
39
So it is feasible but not optimal. Since an probl m with 03 =I 1 will be in£ ibl
we need not fathom that branch, and we are done. The incumbent solution at node
2 represent the maximization of the giv n objective function over th area pecifi d.
4.3 Using CPLEX to Solve Example 2
We can use a similar search technique, with the aid of a computer, to solv the
bouncing ball example introduced ill the previous chapter. This time, instead of
conducting the branchandbound search by hand, we will use the CPLEX libraries.
Thus, this problem will serve as an introduction to using CPLEX and the various
software written to prepare its input and output. For this problem no optimization
occurs, since there is only one feasible solution for each state vector. Vve will be using
the CPLEX MIP solver with a dummy objective function so that it will find the
unique feasible trajectory of the system. A search method using a branchandbound
technique is used to determine which value of b(k) yields a feasible solution for every
k. A flow chart for the analysis is included as hgure 4.4.
A HYSDEL input file which generates an MLD form equivalent to that of Equation
3.29 and 3.29 of Chapter (3) is shown below.
/* bouncing ball
* HYSDEL 2.05
* Kirk Wesselowski
* 31 January 2003
*/
SYSTEM ball_l {
INTERFACE {
/* Description of variables and constants */
STATE {
REAL height [0, 10];
REAL velocity [100, 100J;
}
PARAMETER {
REAL g=9.8;
40
REAL alpha=O.l; /* O=elastic; l=inelastic */
REAL Ts=0.04;
REAL MLD_epsilon=O.Ol;
REAL bounce_t = 0.05;
}
}
IMPLEMENTATION {
AUX {
REAL zl;
REAL z2;
BOOL impact;
}
AD {
impact = height <= bounce_t;
}
DA {
zl {
IF impact THEN height  Ts*velocity
ELSE height + Ts*velocityTs*Ts*g
} ;
z2 {
IF impact THEN (lalpha)*velocity Ts*g
ELSE velocity Ts*g
};
}
CONTINUOUS {
height = zl;
velocity = z2;
}
}
}
Typing the commalld "hysdel ball_i. hys" at the MATLAB command line (the
program can also be run from the DOS prompt) yields the file "ball_i. m" which can
be used to generate the structure "S" which contains all the r levant matrices of the
MLD form, statistic.;s, and error checking results generated by HYSDEL. A campI te
description of the output is given in [28]. Excerpts from the file "balLi.m" are
shown beloVi:
%HYSDEL 2.0.5 (Build: 20020910)
%Copyright (C) 19992002 Fabio D. Torrisi
%%
HYSDEL comes with ABSOLUTELY NO WARRANTY;
%HYSDEL is free software; you can redistribute it and/or
41
!""'""
I. modify it under the terms of the GNU General Public
I. License as published by the Free Software Foundation; either
I. version 2 of the License, or (at your option) any later version.
S.name = 'ball_1';
S.MLDisvalid = 1;
S.MLDstructver = 2;
S.MLDsymtable = 1;
S.MLDrowinfo = 1;
S.Arr zeros(2, 2);
S.Arb zeros(2,0);
S.Abr zeros(O. 2);
S.Abb zeros(O,O);
S.B1rr = zeros(2, 0);
S.B1rb = zeros(2, 0);
S.B1br = zeros(O. 0);
S.B1bb = zeros(O, 0);
S.B2rb
S.B2bb
S.B3rr
°1 °1 ;
] ;
S.B3br
zeros (2, 1);
zeros(O, 1);
[
= zeros(O, 2);
S.Crr zeros(O, 2) ;
S.Crb zeros(O, 0) ;
S.Cbr = zeros(O, 2);
S.Cbb zeros(O. 0) ;
S.D1rr zeros(O, 0) ;
S.D1rb zeros(O, 0) ;
S.D1br zeros(O, 0) ;
S.D1bb zeros(O, 0);
S.D2rb zeros(O. 1) ;
S.D2bb zeros(O. 1) ;
S.D3rr zeros(O, 2);
S.D3br zeros(O, 2);
S.E1 = zeros(10, 0);
S. E2 = [
0.06 ;
9.95 ;
18.01568
17.98432
17.98432
18.01568
190 ;
190 ;
190
42
190 ;
J ;
S.£3 = [
00;
00;
1 0 ;
10;
1 0 ;
10;
o 1 ;
o 1 j
o 1 ;
o 1 ;
] ;
8.£4 =
10;
1 0 ;
1 0.04
1 0.04
1 0.04
1 0.04 ;
o 0.9 ;
o 0.9
o 1 ;
01;
];
8.£5 = [
0.06 ;
10 ;
18.01568
17.98432
0.01568 ;
0.01568 ;
190.392
189.608 ;
0.392 ;
0.392 ;
] ;
S.A = [S.Arr, S.Arb; S.Abr, S.AbbJ;
S.B1 = [S.Blrr, S.Blrb; S.Blbr, S.B1bb];
S.B2 = [S.B2rb; S.B2bb];
S.B3 = [S.B3rr; S.B3brJ;
S.C = [S.Crr, S.Crb; S.Cbr, S.Cbb];
S.D1 = [S.D1rr, S.Dlrb; S.Dlbr, S.DlbbJ;
S.D2 = [S.D2rb; S.D2bbJ;
S.D3 = [S.D3rr; S.D3brJ;
S.nxr = 2;
S.nxb = 0;
S.nx = 2;
S.nur = 0;
S.nub = 0;
S.nu = 0;
S.nyr = 0;
43
S.nyb = 0;
S.ny = 0;
S .nd = 1;
S.nz = 2;
S.ul = inf * ones(S.nu, 1);
S.uu = +inf * ones(S.nu. 1);
S.xl = inf * ones(S.nx, 1);
S.xu = +inf * ones(S.nx, 1);
S.yl = inf * ones(S.ny, 1);
S.yu = +inf * ones(S.ny, 1);
S.dl = inf * ones(S.nd, 1);
S.du = +inf * ones(S.nd, 1);
S.zl = inf * ones(S.nz, 1);
S.zu = +inf * ones(S.nz, 1);
S .xl(1) = 0;
S.xu(1) 10;
S.dl(l) = 0;
S. du(1) 1;
S.xl(2) 100;
S.xu(2) 100;
The structure "S" generated by running "baIL1" at the MATLAB prompt is the
input for the program "mpsoutput.m" written for this thesis. The program generates
the ASCII file "ball_i. mps" which is in the "mps" format. This is an industrystandard
form for writing linear program input file, although it has been adapt d
for MIP files as shown. The program "mpsoutput.m" is complet ly generic with
respect to all types of variables and sy tems, su it should generate valid ".mps" files
from any HYSDEL output. However, each probl m has a specific objectiv function
that must be created. Here we know in advance that this is not an optimizatioll
problem, since the input (gravity) is constant and the system is completely p cified
by the constraints so there is only one feasible solution. We choose a height at
a particular time as the dummy objective function, knowing it will not affect the
solution. Excerpts from the file "ball_l.mps" are shown below in order to describe
its various sections. For a full account of mps format, see [36].
1 he file begins by defining all of the rows of the objective function and constraints.
The first column gives the type of ntry "E" denotes an equality, "L" denotes a
lessthanorequaJto entry, and " " denotes the objective function which has IlO
44
relation operator. The name can be any alphanumeric or p ial
characters; "mpsoutput.m names them ac ording to the typ of ntr ,th tim in th
simulation, and the variable number. Thus, "EROIOOOl' is an equalit constraint at
time k = 1 which refers to the fir t entry at that time. p to k = 99 di crete time
intervals can be tested, with up to 9,999 entries of each time at each tim interval.
ote that the MPS format is a general optimization program form; we are adopting
this notation so we can easily keep track of the trajectory with r p ct to time. Th
first section is shown below, with deleted lines indicated by the ellipsis Lt ••• ' •
NAME HYSDELM
ROWS
E ER010001
E ER010002
L 1E010001
L 1E010002
L 1E010003
L 1E010004
L 1E010005
L 1E010006
L 1E010007
L 1E010008
L 1E010009
L 1£010010
E ER020001
E ER020002
L 1E020001
L 1E020002
L 1E250009
L 1E250010
N OBJECT
The next section names the variables in the program and indic.ates th coefficients
of each variable in each row. Since the matrix has the form Az, with z the vector or
variables, the coefficients are indicated colllmnwise and this is called the "COL M S"
section. The variable name are arbitrary, but here we have CbOSNI to identify them by
the type of variable (XR, XB, ZR, DB for real state, binary state auxiliary real, and
auxiliary binary variables, respectively), the time, and the variable number. Thus,
;(XROIOOOl" refers to a real state variable ("XR"), at time "01" and indicates variable
45
'"""
"0001" of this type at this time. A full account of naming i giv n in th program
listing of 'mpsoutput.m.'· The section marked "I TORG shows th binaryvalu d
variables.
COLUMNS
XR01000l ER010001 001.000
XR01000l lE010001 001.000
XR01000l lE010002 0001.000
XR01000l lE010003 0001.000
XR010001 lE010004 001.000
XR010001 lE010005 0001.000
XR010001 lE010006 001.000
XR010002 ER010002 001. 000
XR010002 lE010003 000.040
XR010002 lE010004 0000.040
XR010002 lE010005 0000.040
XR010002 lE010006 000.040
XR010002 lE010007 000.900
XR010002 lE010008 0000.900
XR010002 lE010009 0001.000
XR010002 lE010010 001.000
ZR010001 ER020001 0001.000
ZR010001 lE010003 001.000
ZR010001 lE010004 0001.000
ZR01000l 1E010005 001.000
ZR010001 lE010006 0001.000
ZR010002 ER020002 0001.000
ZR010002 lE010007 001.000
ZR010002 lE010008 0001.000
ZR010002 lE010009 001.000
ZR010002 lE010010 0001.000
XR020001 ER020001 001.000
XR020001 1E020001 001.000
ZR250002 1E250009 001.000
ZR250002 lE250010 0001.000
MARKOOOO 'MARKER' 'lNTORG'
DB010001 lE01000l 000.060
DB01000l lE010002 0009.950
DB010001 lE010003 0018.016
DB01000l 1E010004 0017.984
DB010001 lE010005 017.984
DB010001 lE010006 018.016
DB010001 lE010007 0190.000
DB01000l lE010008 0190.000
DB010001 lE010009 190.000
DB010001 1E010010 190.000
DB020001 1E020001 000.060
DB020001 1E020002 0009.950
46
DB250001 IE250009 190.000
DB250001 IE250010 190.000
MARK0001 'MARKER' 'INTEND'
The righthand side of each row is sp cified in the eRRS s ction. W pecify
the righthand side set named "N' which say, for example that row EROlOOOl has
a righthand side constant value of 1.92. The first two rows pecify the initial
condition, the rest are inequality constraints. ote that any row with zero on the
righthand side is not included in this section.
RHS
A ER010001 001.920
A ER010002 0000.000
A 1E010001 000.060
A 1E010002 0010.000
A IE010003 0018.016
A IE010004 0017.984
A 1E010005 0000.016
A IE010006 000.016
A 1E010007 0190.392
A 1E010008 0189.608
A 1E010009 0000.392
A 1£010010 000.392
A IE020001 000.060
A 1E020002 0010.000
A 1E250009 0000.392
A 1E250010 000.392
Finally, the bounds on each variable are shown in th "BO DS'" ction. The
and upper bounds can also be specified. The standard default bounc1R OIl any decision
variable x a.re 0 :S x < 00, so if we wish to allow negative values we must specify
notation "FR" denotes a free variable, for example, 00 < XROlOOOl <
the bounds as shown.
. Lower
BOUNDS
FR BND
FR BND
FR BND
FR BND
ENDATA
XR010001
XR010002
ZR250001
ZR250002
47
F"
The file 'balLl.mps" forms the input to the C program te tl. xe which was
written for this thesis by adapting some example file from ILOG that ar fr to u
and modify for CPLEX license holder (such as the MARHES laboratory at Oklahom
State niversity). This C program immediately finds the only £ asibl solution, and
generates two files upon successful completion of the algorithm: "mldmip.lp" how
the form of the program as interpreted by CPLEX, and "MLD_MIP" contains th
objective function value and state vector. The file "mldmip.lp" i in a form that is
easily read by people, as seen in the excerpt below:
\Problem name: mldmip.lp
Minimize
OBJECT: XR200001
Subject To
ER01000l:  XR01000l = 1.92
ER010002:  XR010002 = 0
1E01000l:  XR01000l  0.06 DB01000l <= 0.06
1E010002: XR01000l + 9.95 DB01000l <= 10
1E010003: XR010001  0.04 XR010002  ZR01000l + 18.016 OB010001 <= 18.016
1E010004:  XR010001 + 0.04 XR010002 + ZR01000l + 17.984 DB010001 <= 17.984
1E010005; XR010001 + 0.04 XR010002  ZR010D01  17.984 DB010001 <= 0.016
1E010006:  XR010001  0.04 XR010002 + ZR010001  18.016 DB01000l <= 0.016
1E010007:  0.9 XR010002  ZR010002 + 190 DB01000l <= 190.392
1E010008: 0.9 XR010002 + ZR010002 + 190 08010001 <= 189.608
1E010009: XR010002  ZR010002  190 DB010001 <= 0.392
1E0100l0:  XR010002 + ZR010002  190 DB01000l <= 0.392
ER020001: ZR010001  XR020001 = 0
ER020002: ZR010002  XR020002 = 0
1E020001:  XR020001  0.06 DB020001 <= 0.06
1E020002: XR020001 + 9.95 DB020001 <= 10
1E020003: XR020001  0.04 XR020002  ZR020001 + 18.016 DB020001 <= 18.016
1E020004;  XR020001 + 0.04 XR020002 + ZR020001 + 17.984 08020001 <= 17.984
1E020005: XR020001 + 0.04 XR020002  ZR020001  17.984 OB020001 <= 0.016
1E020006:  XR020001  0.04 XR020002 + ZR020001  18.016 08020001 <= 0.016
1E020007:  0.9 XR020002  ZR020002 + 190 08020001 <= 190.392
1E020008: 0.9 XR020002 + ZR020002 + 190 08020001 <= 189.608
1E020009: XR020002  ZR020002  190 DB020001 <= 0.392
1E020010:  XR020002 + ZR020002  190 OB020001 <= 0.392
ER030001: ZR020001  XR030001 = 0
ER030002: ZR020002  XR030002 = 0
1E030001:  XR030001  0.06 DB030001 <= 0.06
1E030002: XR030001 + 9.95 OB030001 <= 10
lE250009: XR250002  ZR250002  190 DB250001 <= 0.392
lE250010:  XR250002 + ZR250002  190 DB250001 <= 0.392
Bounds
48
•
XR010001 Free
XR010002 Free
ZR250001 Free
ZR250002 Free
o <= OB010001 <= 1
o <= OB020001 <= 1
o <= OB240001 <= 1
o <= OB250001 <= 1
Binaries
OB010001 DB020001 OB030001 OB040001 OB050001 OB060001 OB070001
OB080001 OB090001 OB100001 OBll000l OB120001 DB130001 OB140001
OB150001 DB160001 OB170001 OB180001 OB190001 DB200001 OB210001
OB220001 DB230001 OB240001 DB250001
End
VVe can see from the excerpts above the recursive form of the constraints that define
the trajectory of the system. The "MLD_MIP" file is interpreted as "result_balLi.ro"
and used to generate the plot of the system response shown in Figure 4.3. The
procedure is summarized in the flow chart of Figure 4.4.
4.4 Chapter Summary
Here we have seen how mixedinteger programs are formulated and solved. Many
practical problems are of a ~ize that requires a great deal of computation to solve,
and efficient methods are an area of ongoing I' search. Most methods in lIse today
involve some branchandhound procedure, so it is worth investigating so that we can
understand the CPLEX output presented in this chapter and iIJ Chapter 5, wh re we
will sea more interesting example.
49
•
2
1.8
1.6
1.4
1.2 III
Q;
Q)
E
:i
·wOl
~
0.8
0.6
0.4
0.2
0
0 0.1 0.2 0.3 0.4 0.5 0.6
time, seconds
07 08 0.9
Figure 4.3: Simulation of the MLD boundng ball model solved using the CPLEX libraries.
The coefficient of dissipation is a = 0.1.
50
START
ball1.hys
program HYSDEL at MATLAB prompt
"baIl1.m"
mfile "baIl1.m" at MATLAB prompt
II mfile "mpsoutput.m" at MATLAB prompt
"baIl1.mps"
C program "test1.exe" at the DOS prompt
"MLDMIP" and" mldmip.lp"
FI
Generate plots and analysis using MATLAB.
NISH
Plots and analysis
Figure 4.4: Procedural flow chart for generating solutions of MLD probleml:l using HYSDEL,
MATLAB, and the CPLEX libraries.
51
Chapter 5
Control of a Team of Robots Using
Optimization
The trivial examples of the previous chapters represented certain aspects of the system
behavior in which we are interested, but were too simple to demonstrate how
we can control a true hybrid system. In this chapter, we will provide an application
which show the complete interaction between the finitestate machine and continuous
dynamics for a hybrid system, and demol1strat how we can use an offline
optimization to achieve a desired sy tern behavior. We will control a team of robots
in a simulated game, inspired by the RoboFlag game created by [17] and modeled by
[19], which provides an application for which supervi oryIevel control is of paramount
importance. Our work in this chapter will also help to describe how to use CPLEX
to solve difficult optimization problems, and the CPLEX MIP solver output will be
shown and, to a very limited extent, interpreted.
5.1 The RoboFlag Game
The RoboFlag game proposed by [17] is similar to the familiar children's game "capture
the flag." The purpose of the RoboFlag is to provide a conte t similar to
52
RoboCup Soccer. hut with rules de ign d to promote th d v lopm nt of trategy
and higherlevel control. Current succes ful RoboCup Soc r strategi gen rally put
the emphasis on lower level control, with an individual robot action b iog det rmined
locally. Most RoboCup Soccer teams u e behaviorb ed control strat gie
[18]. The RoboFlag game has more complicated rules and coring so that one can
presume global team strategy will become important. Thus the game is to becom a
platform for developing and testing higher level control schemes. Other l' search fS
have proposed a game of "circular soccer" with similar goals [37].
RoboFlag involves two teams of robots operating on a closed playing field, with
each team initially occupying its own territory of half the field. The object of the
game is for a robot to enter the opponent's flag area, capture their "flag," and retuTll
to their own area of the playing field. A team also earns points for achieving secondary
goals. There is a complicated set of rules for scoring, preventing the opposing team
from scoring and for allowed movements. A full set of rules can be found at [17].
We will model only one aspect of the game, in a simulation similar to that of (19].
The teams are divided into a set of attackers and a set of defenders. The goal of
the attackers is to penetrate th d tense zone of the defend rs (or "score"), while the
defenders move to intercept (or "capt.me" ) the attackers. Th playing field and initial
conditions of a game might look like Figure 5.1, where the d fense zone is reIJI', nted
as a square in order to allow convenient representation uf the zone as a conjunctiolJ
of linear thresholds.
An attacker is captured once a defender has come within a cert.ain distance of it,
thus rendering it illactive. Once an attacker is intercepted, its state latches to th(:_
inactive mode, so the defender is free to move 011 to capture other attackers. The
attacker also becomes inactive once it has scored. The defenders win the game if
they can capture all attackers before anyone of them scores. On the other hand, the
attackers win if anyone of them scores.
53

2 0 2 4 6 8 10
x, meters
8 6 4
1 0 '_'_l~_ __'__ ___>.___'__ __l.___'___ __'___.L__ __J
10
10 I
0 Defender 1
0 Defender 2
8 0 Defender 3
Attackers
6
4
2 0
., 0 D, ~
III 0 E
>
2 0
Defense Zone
4
6
8
Figure 5.1: Playing field for the game, with initial positions for the robots.
54
..
2 0 2 4 6 8 10
x, meters
8 6 4
10 '_'_'__.1..._L_l.__'_'_'__'_'
10
10
0 Defender 1
0 Defender 2
8 0 Defender 3
Attackers
6
4
2 00
l!! D, ~
CD 0
E
>.
2 0
Defense Zone
4
6
8
Figur 5.1: Playing field for the game, with initial positions for the Tobots.
54
We will assume the controller has universal and perfect knowledge of he tate
of the system. As ume the attacker hav no ev ive powerthe proc d dir ctl
towards the goal region at a fixed velocity as long as they are active. A game will be
simulated in discrete time for Tn time int rvals, where Tn is larg enough 0 that all
attacking robots could score if they are not captured.
5.2 Robot Model
Let Nd E IR be the number of robots on the defending team, aIlel let Na E R
be the number of robots on the attacking team. Each defending robot i with i E
{I, 2, ... ,Nd } will have the state vector Xd,i = [Xd,i Xd,i Yd,i Yd,iV, where xd,i E ]R4.
Each attacking robot will have a similar state vector Xa,i = [Xa,i Xa,i Ya,i Ya,i Os,i Oe,i],
where Xa,i E ]R4 X {O, IF. The Boolean variables Os,i = 1 if the ith attacker has
scored, and Os,i = 0 otherwise. Similarly, the Boolean variables 6e,i = 1 if th
ith attacker has been captured by any defender, and Oc,i = 0 otherwise. The state
vector for the entire system of Nd defenders and No attackers is
(5.1)
During simulation, the kinematic input uj(k) for t.lle ,jlh defend r will be controlled
based on the optimization results. The input ui(k) for the 'i''' attacker will be fixed
so that the attacker proceeds to the defense zone as long as it is active, and stops if
it is inactive.
5.2.1 Continuous Dynamics
To model this system, we will use ornnic1.irectional robots modeled in continuous time
as
55
Xi(t)
ih(t)
aXi(t) + (3Ux ,i(t)
aYi(t) + (3Uy,i(t) (5.2)
where (Xi, Yi) represents the Cartesian po ition of the i th robot (attacker or defend r).
The constant a represents a kinematic frictional force on the robot. We use an omnidirectional
model of the robot so that the MIP generated from the MLD form of the
system will be linear. This is an important simplification, since mixedinteger nOlllin
ar programs would be difficult to solve. However, many omnidirectional robot
are used experimentally, and it may be valid for robot used on an actual team.
During the simulation, the i th robot (attacker or defender) will have orne input
series
The input series for the /h defending robot will have Uj determined by solving the
MIP subject to the dynamic constraints presented above, t.he bounds on the tate
and input of the robot, and the logi(:al constraints. In contrast, the input series foJ'
the the i th attacking robot will be constant, so that (ux,i(k), 'uy,i(k)) = (UX ,il Uy,i) for
k E [1, ... , m] unless it is inactive, in which case (ux,i(k), uy,i(k)) = (0,0). The
values of Ux,i and Uy,i will be determined so that the attacking robot travel dire tly
to the center of the scoring region.
5.2.2 Finite State Machine
The finitestate machine of all attacking robot is shown in Figure 5.2. For the i1h
attacker, define an auxiliary binary variable ba,i repre enting the activity tate of the
robot, so that Od,i = 1 if the ith attacker is active, and On"i = aif the ith attacker is
not active. Changes in Os,i and bc,i trigger the transition shown in Figure 5.2.
56
F
~ =1vo· =1 c s
Attacker
Active
8 = 1 a
Attacker
Inactive
8 = 0 a
8 =01\8 =0 c s (8 = 1 v 8 =0) 1\ (8 = 1 v 8 =0) c c s s
Figure 5.2: A finite state machine for an attacker.
5.3 The Robot Model in MLD Form
\Ve will demonstrate the behavior of and interactions between the timedriven dynamics
and eventdriven dynamics shown in Section 5.2 using the concepts listed in
Section 3.2.
5.3.1 AnalogtoDigital Relationships
Capturing
We determine 6e,i E {O, I} based on the r lative posi tions of the attacking and d hlding
robots. In order to maintain the fewest number of linear thr shold relationships,
we will use a square region around each attacker to determine wh n it ha.'3 been captured.
Let 6c,i,j,n be,i,j,l, 6c,i,j,a, and bc,i,j,b indicate whether the requirem nts have beell
satisfied for the ith attacking robot to be captured from the right, the left, abov , or
below, respectively, by the .;th defender. Then for the ith attacker to be captured by
the /11 defender, \ve have
bc,i,j", = 1 ¢::> fe,i,j,r = Xi  Xj  xcap ~ 0
be,i,j,l = 1 ¢::> fe,i,j,l = Xi + Xj  xeap ~ 0
Oc,i,j,a = 1 ¢::> fe,i,j,a = Yi  Yj  JJcap ~ 0
be,i,j b = 1 ¢::> fe,i,j,b = Yi + Yj  Yeap ~ 0
57
(5.4)
where xeap and Yeap ar di tance tha define th aptur tatu us xeap =
0.4 m and Yeap = 0.4 m. If the pia ing fi Id is 20 m wid th n w can d fin
upper and low I' bounds of Me,1' = 19.6 m and me,1' = 20.4 111 for all f, iJ,1" and
Me,l = 19.6 m and mel = 20.4 III for all fe,i,j,l' \ ith a 20 m long pia ing fi Id,
we have Me,a = 19.6 m, me,a = 20.4 m tfe b = 19.6 m, and m ,b = 20.4 m.
Then (5.4) is quivalent to
where f is a small tolerance.
Scoring
Similarl we can d fine 6s,i,)', Os,i,l Os,i,a, and Os,i,b to indicate whether the requirements
(5.5)
fe,i,j,!· < lvl ,1' (1  6e,i,j,1')
.fc,i,j,r > € + (me,1'  € )Oe,i,j,!'
!e,i,j,l < lvIe,l (1  6e,i,j,l)
!e,i,j,l > € + (me,l  €)6e,i,j,l
!c,i,j,u < M ,a (1  6e,~.,.J,a )
!c,i,j,Q > € + (me a  € )6e,i,j,a
fe,i,j,b < McIJ(l  6e i j b)
, " I
!e,i,j,b > € + (me,b  €)6e,i,j,b
for the ith attacking r bot to have scored from the right, the left, above, or below,
respectively. Then
(5 ,i,1' = 1 ¢:} fs,i,1' = Xi  Xdz ~ 0
6s ,i,l = 1 ¢:} fs,i,l = Xi  Xdz ~ 0
6s,i,a = 1 ¢:} fs,i a = Yi  Ydz ~ 0
6s,i,b = 1 ¢:} fs,i,b = Yi  Ydz ~ () (5.6 )
\\ here Xdz and Ydz are distance that define the defense zone, here we use Xdz = 1
In and Ydz = 1 111. We can defin upper and lower bounds of Ms,T = 9 m and
ms,r = 11 m for all !s,i,r, Ms,l = 9 m and ms,l = lJ m for all fs,i,l, Ms,a = D
58
ill and ms,a = 11 m for all fs,i,a, and Ms,b 9 m and m b 11 m for all fs,i,b.
Then (5.6) is equivalent to
zone.
If an attacker become inactive it . hould stop. To model the mode, c1 ctor, we can
5.3.2 DigitaltoAnalog Relationships
(5.8)
Oa,i = 1 t Zi = Uj
Oa.i = 0 t Zi = 0
!s,i,r < Ms,r(l  Os.i.T)
!s,i,r > E + (mS •T  E)Os.i.r
fs,i,l < Ms,l(1  os,i,d
fs.i,l > E + (ms,1  E)Os,i,1
!s,i,a < Ms,a(1  Os,i.a)
!s1i ,a > f + (ms,a  E)Os,i,a
fs.i.b < N/s,b(1  Os,i,b)
f .i.b > t + (ms.b  t)Os,i,b (5.7)
Similar threshold relationships are used to define Od,j,r, Od,j,l, <Sd.j,a, and Od,j,b for the
lh defender so that we can impose the constraint that it cannot enter the defense
define an auxiliary real variable Zi so that
This is equivalent to defining BiCk) = 0 in (2.11) for a particular robot when it is
inactive, and using the appropriate constant B when it is active. Thus, we can reflect
how the state of the FSM affects the switched affine dynamical system. Equation
(5.8) is qui valent to the relationships
mzoa  Zx < 0
N/zoa + Zx < 0
N/z(1 60 )  Zx < ux
mz(1  oa) + Zx < Ux
59
(5.9)
...
where Zx indicate the x component of z for a particular robot and Mz = mz ux '
A similar set of equations applies to the y component of z.
5.3.3 Logical Relationships
Capturing
We know that the i th attacker has been captured by the ph defender when
be i] = 1 B be i ]' r = 1 /\ be i ] I = 1 !\ Oe;]' a = 1 /\ Oe i ]' b = 1 " " , " t ,", t· " ,
(5.10)
then
Similarly, V\le know that the i th attacker has 1;cored when
0 ,. + 0 ' ' C,t,],r C,1,]
...
(5.11)
0e,"t,]I, + 0e,l',],
bc,i,j,a + Oc,i,j <
 Oc,i,j,b + Oe,i,j <
be,i,j,r + Oe,i,j,l + Oc,i,j,a + Oc,i,j,b  Oe,i,j < 3
< 0
< 0o
o
Scoring
(5s ,i = 1 B Os,i,r = 1/\ Os,i,l = 1 /\ Os,i,a = J /\ 0s,i,& = 1 (5.12)
then
<S5,,t,T +0S,t',l+OS,'t,a +O8,1,bO8,'1 < 3
fJs i r + fJs i < 0 " ,
Os,i,l + Os,i < 0
Os,i,a + Os,i < 0
Os,i,& + Os,i < 0 (5.13)
60
where Oe,i is determined by the relationships in Section 5.3.5. The inactive state is
terminal because of the FSM (or Automaton) definition.
......
Activity
For the ith attacker, the activity is controlled by
which is equivalent to
Os, i+Oa,i < 1
oe,i + oa,i < 1
Os,i + oe,i + oa,i > 1
5.3.4 Dynamics
When we discretize (5.2) and model the MS as in Section 5.3.2 with Q
(3 = 1 at Ts = 1 sec for the i th robot, we get.
[
1 0.1 lJ
o 4.54· 105 0
o () 1
000
[
lJ.09 0 ]
.
0.10 0 .(k) + 0 0.09 Zt
o 0.10
(5.14)
(5.15)
10 and
(5.16)
..
...
The general form of the discretized statespace representation of the system is
x(k+l)
y(k)
Ax(k) + B1u(k) + B2o(k) + B3 z(k)
Cx(k) + D1u(k) + D2o(k) + D3z(k) (5.17)
where u(k) = [u;' ut]T is a vector of inputs, and y(k) = [y;' ytV is a vector of
outputs. By using the appropriate matrices in (5.17), we can realize the behavior of
61
the SAS as described in Section 5.3.2.
5.3.5 Autonlata
The automaton section defines the behavior of the finite stat rna hine d crib d in
Section 2.2.1. The i th attacker can be captured by any d fender so th Boolean tate
variable .5c,i is controlled by
individual latches on the captured state and scoring state, instead of latching tIl('
Figure 5.2 shows that the inactive condition is a terminal state, thus once the FSM
enters the terminal state it "latches" to it. In order to facilitate debugging, we'll put
Nd
.5c,i + L 6c"ij < 0
j=l
<5 ,.5 1 > 0 C,t C,'l.,
<5 . .5 . 2 > 0 C,'l c,t,
15e,l. .5C,t.,Nd > 0 (5.18)
activity variahle. For the captured state, we lllllst have <5c•i (k) t .5 ,i(k + 1) at tilll('
interval k, or
(5.19)
where we show the time interval k tbat had heretofore been suppressed. The automaton
equation that achieves this, while controlling the incoming edge defined in terms
of the variables of (5.11) is
(5.20)
Similarly for the scoring state, we must have .5s,i(k) t os,i(k + 1) at time interval
k, or
62
!"""
6s,i(k)  6s,i(k + 1) :::; 0 't/i E {1, 2 ... ) a}
5.3.6 Constraints on the State of the Defender
(5.21)
Since no defender can enter the defense zone, we have to constrain Xd,j for the jLh
defender. Recall the definitions of 6d,j,Tl 6d,j,l, 6d,j,a, and 6d,j,b from Section 5.3.1; then
the literal statement
must be true, which is equivalent to
which can be implemented using the "MUST" section of HYSOEL.
5.4
6dl]')r I\ tJdIJ'~ll\ tJd,)',a I\ tJd,J'b,
(6d' v6i 'I) V (rd , VtJd 'b) ,J,r ( ,J, ,J,a ,J,
Game Simulation Results
(5.22)
(5.23)
In this section we describe how to solve the problem by building the PLEX iUjJnt.
file in MPS format and running th CIanguage command file. B cause the problems
to be considered are mu 'h more substantial than the cxamples prcviously shown, it
will provide a good opportunity to discuss how to interpret the CPLEX output.
The HYSOEL input file "Robo3def8att.hys" uses the relationships c1escrilJed above
in the "LOGIC,' "AD," 'lOA," "CO TIN OUS," "AUTOMATON," and "MUST"
sections of the input file. The MPSformat input file generated by the MATLAI3
function is similar to the examples shown in Chapter 4. The objective function to be
minimized is
N"
V(k) = L tJs.i(m)
i=]
63
(.5,24)
MINIMIZE
OBJECT
A
BND
A typical CPLEX session output u ing the default parameter ttings i hOWD blow.
It shows that even the default mode us a impl 'PI'  olv stage that eliminat
redundant rows.
Selected objective sense:
Selected objective name:
Selected RHS name:
Selected bound name:
Tried aggregator 3 times.
MIP Presolve eliminated 1984 rows and 718 columns.
MIP Presolve modified 76 coefficients.
Aggregator did 982 substitutions.
Reduced MIP has 9174 rows, 3694 columns, and 24846 nonzeros.
Presolve time = 0.22 sec.
Clique table members: 2664
MIP emphasis: balance optimality and feasibility
Root relaxation solution time = 1.02 sec.
Nodes Cuts/
Node Left Objective lInf Best Integer Best Node ItCnt Gap
0 a 0.0000 294 0.0000 1958
0.0000 427 Cuts: 603 3041
* 0+ a 0 2.0000 0.0000 3041 100.00%
* 10+ 9 0 1.0000 0.0000 4039 100.00%
100 18 infeasible 1.0000 0.0000 6031 100.00%
* 130+ 1 0 0.0000 0.0000 6448 0.00%
Flow cuts applied: 72
Gomory fractional cuts applied: 80
Solution status 101.
ans : <
Solution status = 101>
ans : <Solution value 0.000000
>
Objective value a
This solution was generated in 27 seconds (actual time) using a PC with a Pentium
4 processor. It represents an optimum solution according to the objective hlllction
(5.24). Vve know it is in fact globally optimum, since we did not give CPLEX any
criteria other than the final score as a measure of performance, and the lowest possible
score IS zero.
The output shows the following information:
64
F
Node  The solution status is printed for very integer ~ ible olution indicated
by a *, and for every lOoth node of the optimization tree, as shown above.
Nodes Left  This column indicate the number of nodes not xamin d.
Objective  The objective function of the relax d solution is hown, unles th r
is no feasible solution even for the relaxed problem.
Iln!  An integer feasible solution has no integer infeasibiliti ; otherwi e, the
number of infeasible binary variables (that do not have a value of 0 or 1) is ShOWll.
Best Integer  Shows the objective fun bon value for the current incumbent solution.
Cuts/Best Node  The current objective function value is compared to the be t
objective function value for the relaxed problem, subject to any applied cuts. Thus,
the Cuts/Best Node value starts at the objective function value of the original relaxed
optimization prohlem (the root node), here 0, and may hecome worse as cuts are
applied. The total number of each type of cut applied is listed at the end of the
output.
!tent  Shows the number of iteratiolls so far, including all iterations at each
node.
Gap  Shows the difference between the urrent solution and the uts/Best ode
column.
The resulting simulation is shown in Figure 5.3. Note that we have demon trated
the behavior of a true hybrid system while simulating the cooperative control team
of robots. Each defender avoids the defense zone &'5 required, which we achieve using
t.he auxiliary binary variables and EG as shown above. Each attacker stops when it
is captured, 0 the FS tI affects the SAS through the MS as required. The defender
kinematic inputs are limited to I U.,i I:::;: S m/s as shown in Figure 5.4. So even though
our model is simple, we have demonstrated the ability to conform to constraints on
the continuous and Boolean states of the system and find an optimum solution.
65
.1
..
10
8
6
4
2
~
.l!!
III 0
E
>.
2
4
6
8
10
10 8 6 4 2 0 2 4 6 8 10
x, meters
Figure 5.3: A trial game with 3 defenders and 8 attackers. Bvery ;j,tlacker (in reel) IS
intercepted. by a defender (in blue, green, or purple) before scoring.
66
..
i sr==!:::::;;::::::;~::::;rr"'';~r=='=H
'5
Q.
Eg 0
e
c:
lsL.:~_____J~::::::;~::::j==:::::;:~=:i:==::::::=____'_~~_..::..._.J
o
i SL=='==.,.,.:==:r.::r:;==::;]
'5
.~
~ 0
e
~
~ 5 L_~_~==::::=:::..._ ___..:::__..L..:._'__ ____L.:::._1::~.J
o
Figure 5.4: Kinematic inputs for a trial game with 3 defenders and 8 atttl(.:k~rs.
Of course, it might be advantageom3 in a real game to refine our notion of optimal
performance. For example, we may wish to minimize the score wlJile using the lowest
po sible amount of input energv, which we ma.y achieve by rninimizing the sum of
kinematic inputs to tIL defender' over the course of the simula.tion. W(~ 'an a.chiev
thi by minimizing the objective function
" Nd
V(k) = I>~3,i(m) +WuE IIUjih
i=l j=!
(5.25)
where Wu is a weight for the input terms and Uj is defined in (5.3). The nonlinear
onenorm fUllcti.on is achieved by setting the elements of an auxiliary real variable Zu
equal to the absolute value of the element of u based on an auxiliary binary variable,
using a method similar to that shown above.
We can interrupt the branchaDdbound at a preset maximum tim if desired, at
which point CPLEX returns the incumbent solution as the optimum. For this trial
we limited CPLEX to 20 minutes of actual computation time, and at 20 minutes the
67
..
2 0 2 4 6 8 10
x, meters
4
. 0
'. ' ... O.O8g000o0~ oo .
0.' ... 8·· ...
.ooooOOQ ..
10 '_'_'_L__'_L_'_'__J_~_ _'
10 8
10
8
6
4
2
~'"
41 0 E
>.
2
4
6
8
The solution does appear to corne close to successfully minimizing input ClI rgy for
Figure 5.5: A trial game with inputs minimized. The defenders are successful in preventing
scoring, but use less energy than in the trial shown in Figure 5.3.
solution shown in fi'igur 5.5 a.nd 5.6 was rcturned. Ther [or , adding the r fJuirC'm
nt that input energy bc minimized significantly incr a' d th omputational I act.
example the integTation of total kinematic input over the cour e of the simlliati Tl is
498.8 m/s for Figure 5.4, but only 150.4 m/s for Fignre 5.6. ote that ill our problcm
formulation, the capture zone if; a box of width and height 0.8 TIl, but the captured
attacker does not stop until the time interval after the one in whicL the criteria for
the state transition were satisfied, because of the dynamic FS model formulation.
Thi makes the capture zone appear larger than its actual ::,'ize in Figure 5.5.
68
~.. 5
;;
a.
S
lS 0
~..
~5
0 2 6 8 10 12 14 16 18 20
time
~'" 5
S a.
.S
.8 0
2
c:
L s '" 0 2 6 8 10 12 14 16 18 20
time
co
~ 5
Sa.
.S
.8 0
2..
~ 5
a
0 2 " 6 8 10 12 1<4 16 18 20
~me
Figure 5.6: Kinematic inputs for the trial game witb inputs minimized.
5.5 Chapter Summary
Even for the simple system shown here, CPLEX solution times on a personal eomput'r
an be slow. It is important to formulat the pr bi m in an effiei nt wa and to
ontrol the CPLEX solver in order to achi v useful solutions in a timel.y HI.ann 1'.
We have demonstrated how we can control of a hybrid system in all offlin imula i TI,
and pointed out some of the interesting salient featnT<*l. Our II xt topi will b tIll!
imulation of a recursive, online control system, again using optimization t chlliqlle~.
69
Chapter 6
Model Predictive Control of Robot
Formations
The pre\ ious chapters have described control of robot teams by first describing a
system ha.ving eventdriven dynamics and continuous dynamics, then solving the offline
optimization problem to maximize a given objective function. Here we consider
model predictive control (MPC) for the online feedback control of coordinated robot
formations. Because many mobile robots arc modeled as nonlinear systems with nOrtholonomic
constraints, the inherent ability of MPC to handle nonlinear constrained
sy tern might make it a good technique for Dlultirobot formation control tasks such
as trajectory tracking. In this chapter, a simple MPC algorithm for control of three
unicyclemodel robots tracking a formation is simulated. Th formulation us s a
terminal cost to create stability.
6.1 Predictive Control Methods
Passive state feedback laws have been well explor d for control of nonholonomic robot
systems[38], but recently some research has applied modelpredictive cont1'Ol (MPC)
to control of formations of robots [20]. Broadly speaking, MPC is class of algorithms
70
...
'"""
that rely on an optimization to determine th be t alu for con rol inpu in rd r to
achieve specified goals. This optimization require a model for the t m (or plant)
being controlled, and actually generates a pr diction for the futur r pon e of th
system. This is in contrast to state feedback eontrollers that do not consid r a model of
the system. In such controllers, the ouly "feedforward" component is in the d rivative
term of a PID controller, which is in a sense an empirical feedforward compon nt.
However, the prediction inherent to MPC gives the method at dforward component
that will allow the system to respond to known and anticipated input di turbanc s
or changes in the reference trajectory. At the same time, the sy tern respond to
measurements of the actual state of the system, allowing a feedback re pon e to
input noise. This combination of feedforward and feedback components may allow
MPC to outperform "passive" state feedback controllers in some applications.
State feedback PID controllers are applicable to linear systems, although extensions
such as gain scheduling may allow the basic concept to be applied for nonlinear
systems. However, MPC can be applir.d to nonlin ar plants subject to hard and soft
constraints. The ability of MPC to illrorporate hard constraints (that cannot be violated)
and soft con traint (that can he violated but with som p nalty) is another
distinct advantage of MPC. Plants with state fe dback PID controllers r pond to
Gaussian input noise with some Gaussian output, hopefully with an acc ptably small
variance. Because disturbances on either side of the refer nce trajectory are handled
with symmetrical corrective action, the plant must be operated at some distance from
any hard constraint. In contrast, it is sometimes possible to operate plants controlled
by IPC doser to hard constraints, because the optimization that is the heart of the
MPC algorithm can be a nonlinear algorithm that can find solutions close to hard
constraints. In other words, it does not need to respond symmetrically to noise on
different sides of the mean value, and Gaussian input noise can have a nonGaussian
output [39].
71
...
These are some of the advantage of MPC that have au d it t b om increasingly
popular for control of complicated terns uch in the p troch mi al
industry [40]. However. it. doe have disadvantages. In it traditional u for process
control, the primary disadvantage is often considered to b the n d for a good
model of the plant. In application to robotics, the foremost disadvantag may be th
computational cost, which is often negligible for slowmoving sy tern in th pro e
industry. Because nonlinear optimization is an iterative earch algorithm, it can b
a slow process for large plant models which can ea ily have thousands of con traint .
The optimization of systems \vith nonlinear objective functions and / or nonlinear
constraints is generally not guaranteed to find the globally optimum solution. On
the other hand, delays in control system calculations can be detrimental to system
performance. Since the optimization routine achieves progressively better solutions,
an important area of research is in the interaction between the realtime embedded
system controller and the realtime plant [41]. However, the potentially high computational
cost of MPC has often limited its practical applications to those with slower
time constants.
Although the MPC formulation seems quite intuitive, using a seemingly commonsense
formulation without analyzing stability can lead to divergent re pons s. It
may not be difficult to prove tability for infinitehorizon control laws, or for MPC
algorithms that require a zerostate constraint at the end of the horizon at each tep.
These approaches often do not lead to tractable problems. With shorter horizons that
can be handled computationally, a stabilizing solution may not exist. Much of the
recent research in MPC has centered in analyzing stability using Lyapunov analysis
centered on appropriate terminal costs or constraints. See, for example, Mayne, et.
oJ for a good survey of recent MPC research [42]. Certainly, a stability analysis
should consider what constitutes a meaningful definition of stability, and practical
methods of ensuring it.
72
•t•
)
)
With improved computers and calculation m thod it rna b po ible to u e
MPC for control of ystems with generally fast time con tant uch as in mobil
robotics. This could be an important area of re earch because the advantage of MPC,
particularly the ability to handle constraint and nonlin ar mod I ,ar well suited
t.o solving control problems for formation of robots. Recently, ome imulations
have been developed to explore the use of MPC for trajector tracking of formations
of robots. We will demonstrate an MPC algorithm for the trajectory tracking of a
threerobot formation in an equilateral triangle formation, similar to the method of
[20].
6.2 Model Predictive Control Formulation
Consider a dynamic model for a system to be controlled in the form x(k + 1) =
f(x(k), uo(k), .6.u(k)) so that f(·) is time invariant. and. nonlinear. Assume that
the input vector is possibly constrained with some upper and lower bounds so that
Umin :S u(k) :S U max' We could drive output of a system from some initial
condition y(k = 0) = Yo towards ree renee trajectory {r(O), r(1), ... r(kmax )}) where
r(k) E jRH is some desired value of the state vector at time k, by applying somp.
appropriate sequence {u(O),u(I), ... ,u(kmax I)}. L t u(k) = uf(k) + .6.u(k),
where uf(k) the input required to maintain zero error in the system trajectory, awl
.6.u(k) is applied t.o correct any errors. We want to minimize the cliff ren 'e b tween
y( k) and r( k) as time progresses. In order to have a realistic and economical control
law, we also want to minimize the amount of input energy required to achieve this
goal. Here, let any control input series that stabilizes the errur between rand y
to the origin in Hp steps or fewer while not violating any constraints be denoted
u(k) = (.6.u(k), .6.u(k + 1), ... ).6.u(k + Hp )), and let an optimum control input.
series that stabilizes the system in the least possible time and/or with the least
73·
..
possible effort (while not violating any con traints) be d noted iLO(k).
The objective function V(k) = g(x(k), iL(k)), will be generall proportional to
control effort and error of the predicted trajectory from the d ired traj tory. Th
objective function has the general form
Hp
V(k) = L p(x(k + m), 6.u(k  1 + m))
m=l
(6.1)
where pC) is an incremental cost. An MPC formulation is solved u ing an optimization
routine run over a finite prediction horizon, B p , and appli simulation inputs over a
input horizon, B u , where Bu ::; B p . The optimum objective function value associated
with ij,°(k) will be denoted VQ(k). An optimization routine varies th control input
over some control horizon Hu , where B u ::; B p •
A generic MPC response for a singleoutput system is shown in Figure 6.1, where at
time ko (zero on the abscissa) the output y(ko) is in some error from the desired s(ko).
This error could be a result of some system disturbance, error in state estimation,
etc. We wish to determin orne input series 6.u(k) = uc(k) for ko ::; k ::; Hu
which is to be added to the original free input uf(k) to return the sy tern to th~
desired reference trajectory. For any time intervals BtL ::; k ::; Hp  I, the input
vector is typically fixed at the value u(Hu ) as shown in Figure 6.1. The d sired
trajectory r(k), ko ::; k ::; B p [or correcting the state of the system is often an
asymptotic approach to the setpoint trajectory s(k) (howev r, for our robot system,
they coincide). Because the system model in general may not allow this firstord r
response as shown, the predicted trajectory Yp( k) may not coincide with the r ference
trajectory r(k). Although Figure 6.1 shows a constant s(k) and uf(k), this is not
necessary.
TllP optimization routine will find the uc(k) to minimize the error between Yp(k)
and 1'(k) subject to the appropriate constraints. While the state vector evolution is
C'onstrained by the dynamics, the input vector is free (subject to the given upper and
74
,.
)1
,J.'
~  
ooeoQGOe
0.8
s(k)
00 0 \
y(k)
~06 yp(ko+Hp I~)
r(k) ~ 0
Yp(kO+Hu I koJ
'R0.6 " ~OO
:;
/0
0
0 0.4 / /0 0 "Yp (klkOl
0.2 /00
\Y(kol H H
0 u p
10 5 0 5 10 15 20 25 30 35 40
time, k, as difference from current time step ~
ue(kolkol
...0
0
0
0
o U (klkol
o ~ e
00
00 ue(ko+Hu I kol e(kO+Hp1 I kol u,(k) 00 ~ooooooooooooj J 00
Hu Hp
0.8
:; 0.6
0.
.£
0.4
0.2
o
10 5 o 5 10 15 20 25
lime, k, as difference from currenltime slep ko
30 35 40
Figure 6.1: General form of an MPC response. The reference trajectory T'(k) is planned
to return the system to the setpoint trajectory s(k), but the predicted response Yp(k) may
differ from T(k) because of latencies or const.raintR.
lower bound) to be adju ted by the optimization algorithm so as to best. achieve
Ollr objective given hy (7.4). At each discrete time interval k, an optimization is
performed, but only the first of t.he sequence is actually implemented. Instead of
irnplementing the remainder of the optimal sequence {u(k + 1), ... ,u(Hp )}, w run
another optimization routine at the next time interval, using the actual state at
that time interval. This allows the system to respond to any timevarying or fix d
disturbances. The value of k at Hp then increases at each time, so it is a receding
horizon. Note that we have implicitly assumed that u(k) is available at time k, so
that calculation times are n gligible.
75
1
., ,"
)
to:
6.3 MPC for a ThreeRobot Formation
Although many formulations of t,IIPC have been develop d [43], w will u a stat 
space representation since this is a common way to present robot model. Consid r
a state space model for robot i at time k such as
Ii (Xi (k), Ui ( k) l k)
9i(Xi(k), ui(k) k) (6.2)
(6.3)
where xi(k) E X is a state vector, ui(k) E U is an input vector, Yi(k) E lRTn is an
output vector, and the functions I and 9 are Coo on X and U respectively. \Ve will
consider unicyclemodel robots with the continuoustime kinematic model
r
:~(t)] rc?sB(t) 0] [vet) ],
~(t) = sm (J(t) 0 ' (t)
B(t) 0 1 w
where x(t) E IR and yet) E lR are the Cartesian coordinates of the robot, and I" 
B(t) E IR is the angle of the robot with respect to the positive x axis such that.
7[" < (j(t) :S 7[". We ,;\,rill use a discretetime representation of (6.3). At any discrete
time k the state of the robot is
xa(k) = [x(k) y(k) (} :i:(k) iJ(k) O(k)]T (6.4)
),)
..
The model inputs are taken to be the linear velocity v(k) and angular velocity w(k),
and we do not consider any dynamics required to convert an actual motor torque
to these velocities. The unicycl model in (6.3) will represent the kinematics of all
robots. At any discrete time k the desired state of a robot is
(6.5)
VVe would like to drive the actual state of the robot X a(k) to xd(k) in some stable
manner, so that the error vector x(k) = [xd(k)  xa{k)] + 0 a') time increases.
76
..
The state of the sy tem x(k) should be an element of a conv x, closed t X c ]R6
containing the origin. The state is can rolled by the input v ctor u(k) = uf(k) +
.6.u(k) = [Vj Wj]T, where uf = [Vo wo]T. The function Vo and Wo repr nt the
velocities necessary to maintain zero error (i.e. de ired formation), and are bounded
if the kinematic inputs of the leader robot are bounded. Thus,.6.u = 0 when
x = O. Each .6.u(k) is an element of a convex, compact et U C ]R2 containing the
origin. Vve will define U = {.6.u(k) I Umin :::; .6.u(k) :::; u max } so that th kinematic
inputs are limited to some reasonable magnitudes. For the work in this chapt r, the
plant will consist of three kinematic models, one for each robot, running concurrently.
Consider Hu = 5 and Hp = 6, and let the vector of decision variables be z(k),
where
...
,
I
xs(k + 1)
z(k) = us(k + 1) (6.6)
...
.i". ,
l
I...
xs(k + 5)
us(k + 5)
x s(k+6)
where xs(k) = [xf(k) xHk) xnk)]T is a vector of stat variables for th thre
robots, and us(k) = [u{(k) uf(k) uf(k)F' is a vector of input variables. The valu
of xs(k) will be fixed at xo, defining the initial condition for the given iteration at
time k. The desired pose of the robot formation consists of static r lations between
the robots, and a dynamic movement of the formation along the desired trajectory
r(k) = [rx(k) ry(k)jT. Therefore, if the plant were operating with no deviation
from the desired evolution, the terms defining the interrobot separations and angles
would not change, and the whole formation would move as a unit along the desir d
trajectory. Let the desired distance between the center of each robot and the point. 011
the reference trajectory be a constant d, then with the robots forming an equilateral
77
.
j:,)
,.
\
Reference Trajectory R
\ "
\
\
\
\
31/2 d
d
/
"
I (R (k),R (k)) /
• ll; Y /
I /
I / 3
112
d
Id
\
:J ........ ...
Figure 6.2: A 3robot formation tracking a reference trajectory in an equilateral triangle
formation. The arrow shows the robot orientation.
triangle the interrobot separation will be V3d as shown in Figure 6.2.
The development below follows that of !20], except in some details. The objective
function will have the forIll
.. ... ,.
"II
II
,~.
".
(G.7)
r:
'I
'1
.'I.
where the objective function cost terms Ep and Es penalize any deviation from some
desired position and separation, respectively, and Ef(k) is a cost associated with a
deviation from a desired formation. Here we will use
Ep(k)
fj 3
2:I: Wp(i)
i=] j=]
(J(Xj(k + i)  xs(k + i))2 + (Yj(k + i)  Ys(k + i))2 + ~  d)2
6 3
2: I: Ws(i)
i=] ),1=1
j#/
78
P'"""
where the prediction horizon Hp = 6 Wp(i) and WS(i) ar w igh at time i for th
position and separation terms. respectively, and ~ is a small arbitrary constant. The
desired state Xd is implicit in (6.8).
In addition to the proper separations, the optimum traj ctory for the controlled
system should use input power that is minimized according to some weight ltVu(i), 0
we consider the terms
5 :3
Eu(k) = L L WU(i)u(k + i)
i=O j=l
(6.9)
The formation achieved with the objective function in (6.8) would not have a stable
orientation around the reference trajectory. Therefore, we must consider an objective
function term that penalizes a formation deviation from the desired shape and
orientation. To achieve this we will add a terminal cost to our objective function.
The formation orientation defines a desired equilibrium formation consisting of a
Sf't of desired permissible states for the robots when they are following the trajectory.
Since any of the three robots can take any of the t.hree positions, the equiliurium
formation is precise but not unique. The formation is achieved according to (6.8)
when Ep and Es go to zero. This formation is not stable, jn that there is a continuum
of f asible states (imagine the equilateral triangle formation shape rotating about
r(k)).
To implement the formation orientation we must use heuristic rules to define
which robot should adopt each position in the formation. In the MATLAB code used
to implement the MPC algorithm, the program determines the norms of errors in
the Cartesian coordinates of the robot from the desired positions for every possible
combination of assignments, then makes the assignments that minimize the sum of
these norms for the three rohots. An error vector of the form
'.
'.
I
h'
'.;::
III
(6.10)
79
is used, and the terminal co t terms whi h d fin the formation orientation hav the
form
3 6
Ef(k) = L L x;(i I k)Wfxj(i I k)
j=1 i=2
where Xd,j and Yd,j are the desired Cartesian coordinates for robot .i and Wf i a
constant diagonal weighting matrix. The MPC formulation is
mznzmzze
with respect to
subject to
x(k)
xl(k+1)
x2(k + 1)
x3(k + 1)
x](k+6)
x2(k + 6)
x3(k + 6)
u(k) <
u(k) <
6.4 Simulation Results
V(k)
u(k)
Xu
11 (xl(k), u\ (k))
.!2(x2(k), u2(k))
f3(x3(k), u3(k))
11(Xl(k + 5), ul(k + 5))
12(x2(k + 5), u2(k + 5))
h(X3(k + 5), u3(k + 5))
u max
Umin (6.11)
...
.:
r~
\I
I
A simulation program was developed using MATLAB. Although MATLAB has a commercia]
MPC toolbox, the algorithm did not use it. This system is not complicated
and the algorithm was generated specifically for this project by taking advantage
of MATLAB's "fmincon" nonlinear optimization algorithm from the Optimization
toolbox.
A typical simulation is shown in Figure 6.3. Tote that the relative positions of the
red, blue, and green robots depend on their initial (;onditions. Thus, any of the three
80
W'C of three nonholonomic robots
'.2
!
.!l
~ G G ,:;
G G
0.5 '.5 2 2.5 3 3..5
lC, meters
t'!
~
E ,;, G G
0.6 G G
0.4
0 0.5 1.5 2 2.5 3 3.5 ..
... meters
1.2
G
l! G !! 0.8 ~ G G G ,:; 0.6 G G G
0.4
0 0.5 1.5 2 2.5 3 3.5 4
... metern
ft'igure 6.3: Evolution of the 3 robot fonnatioD following a trajectory from different init.ial
conditions.
can end up in the bow port, and starboard positions (we use the nautical tenns to
emphasize the lack of I ader and followers).
Our simulation shows that if the terminal cost is neglected by setting tIle w i hts to .,
zero, the formation 1,a.'3 no stable orientation as shown in Figure 6.4. he re pon'
is slower, but the system eventually sta.bilizes to a certain set of 10 'ations arouTlCl
r(k) (off the plot). This demonstrates a certain level set of X, but not an ace ptabk
stability.
The formulation with the terminal cost does a.ppear to eventually stabilize from
allY la.rge initial errors tested. Figure 6.5 shows tht responsc tra,(~king a line trajectory
from the initial conditions :1:1 [0 1  7rjT, Xz [3 1.5 7r/2rf', and X3
[0 0 0]1'. The upper plot shows the response when the robot i prevented from
driving backwards. This restriction is arbitrary (and the simulation i. not realistic),
81
2.5
2
.I!.! 1.5 " " . 1
ii • we w u
E 0.5
>.
0 "
<l.5
1
0 2 6 8 10 12
" melers
2.5
2
1.5 " .~.  ~ •
0.5
>.
0
<l.5 )
1 .
0 2 6 8 10 12
" meters
Figure 6.5: Evolution of the 3 robot formation from larger initial errors. The upper plot,
shows a system where the robots cannot have a negative input linear velocity; no sud.
restriction is used for the lower plot.
Our simulation is not run in real time, but important area for future research could
be to implem lit a realtime simulation so that th cost of PU tim' is in 'orporatcd
into the evaluation of the algorithm. ft rna be possibl to int Tmpt tllC algoritluTl
based on some preset critcl;a, balancing the need for input eommands at re lIar
intervals with the de ire for optimalityI41].
6.5 Chapter Summary
We have seen that certain requirements must be met for the system to exhibit good
stabilityin particular, the sy tern is not stable if we do not applv the termim),]
cost show·n. Because we want to understand different approach to rnodel predictive
control, we have also explored another approach to the problem u ing a dualmode
model predictive control scheme. This is tIle subject of Chapter 7.
83
...
...
Chapter 7
DualMode Model Predictive
Control of Robot Formations
In this chapter, we present another method of establishing stability for MPC systerns;
that is, using a dualmode model predictive control algorithm with a terminal
constraint set. We develop a dualmode controller for teams of nonholonomic robot
in leaderfollower configuration and contrast it with the use of an inputoutput fe dback
linearized controller alone. We discuss the advantages and disadvantages of
dualrnode control versus other predictive controllers ill thi. application, and show
how the concept may be extended for control graphs of larg r arbitrary formations.
Conditions for stability are discussed and analyzed.
7.1 LeaderFollower Algorithms
'vVe will compare and contrast two methods for formation control. First, th(~ leaderfollower
approaches involve a team of robots following a real leader robot or a virtual
robot. Leaderfollower approaches may be susceptible to failure of the leader or leaders,
or may suffer from a lack of feedback to the leader if the follow rs are impeded.
Second, decentralized methods, often based on optimization problem , may be more
84
)
,
."

robust (with regard to individual robot performan ) and mor fl xibl in m ting
performance requirements, but may be computationally expen i. de ntraliz d
approach is taken in [14]' where the lack of robu tne of the laderfollow r m thods
is overcome by using clever elementary formation moves. Formation b havior is
driven in a series of formation moves that are assum d to be feasible and re ult in no
collisions.
The coordination of upperlevel or supervisory control with low rIevel or local
control has been an important topic for formation control [44] [16]. In [44], control
graphs are defined that allow the separation of control tasks into group, formation,
and local control levels. A heuristic algorithm is proposed in lieu of a mixedinteger
program to to optimize the formation control level. In [16] stability is shown for
control of fullyactuated (holonomic) vehicles, using a method with virtual leaders
and potentials.
One approach for control of formations of robots is to use an inputoutput feedback
linearization control law applied t.o a "follower" robot at some sp cifi d relation hip to
a "leader" robot [11]. sing this technique, complex formations can be achieved with
defined relations for chains of robots. Each robot in the formation rna b following
one robot and leading one or more other rohots. Stability of formation control using
this method can be proven under appropriate conditiom; [45]. Oth r researcher hav
provided somewhat similar algorithms using virtual leaders [13].
In contrast to Chapter 6, here we will use terminal constraints to en ure stability
of the MPC algorithm.. From the terminal constraint region, we use a stabilizing
local controller in the manner of [23] and [24]' resulting in a dualmode MPC algorithm.
This may have the advantage of conceptual simplicity, hut may suffer the
disadvantage of infeasible optimization problems outside a necessarily limited region
of convergence for the algorithm. We contrast the dualmode MPC controller with
the use of an inputoutput feedback linearization control law alone. In the dualmode
85
po
MPC algorithm, we can incorporate a colli ion co t term to r duce or if tiv 1
eliminate the possibility of interrobot colli ions. Eith r control m thod may form
the lowerlevel control for coordinated formation control sy terns.
This chapter is organized as foHows. In Section 7.2 we d scribe the probl m formulation
and how it can be solved u ing inputoutput feedback linearization controller
and dualmode MPC algorithms, and in Section 7.4 we how simulations contrasting
the performance of these two methods. In Section 7.5 we describe how these methods
could be used as the lowlevel control with an overarching graph theory approach to
formation control. Finally, in Section 7.6 we provide our concluding remarks.
7.2 Problem Formulation
vVe will again consider the planar motion of nonholonomie unicyclemodel robots
whose kinematics are determined by (6.3).
7.2.1 Inputoutput feedback linearization
We will use an inputoutput feedback linearization controller as d(~vcloped by [11],
where full details can be found. For cOlltrol of two rohots as shown in Figure 7.1,
denote the leader robot as ~ and the follower robot as Rj . Since ware using a
timeinvariant contimlOllS feedback law, we control the point Pj which lies a distance
c1 normal to the center of the axle of the robot. Consider that the leader moves
independently given some arbitrary input Ui = [Vi WijT. The follower Rj is controlled
so that its trajectory approaches a position rt t some desired distance ifj and desired
orientation 'l/Jt relative to the leader. The control law determining Uj = [Vj Wj]T
based on the position of ~ has been denoted a Separation Bearing Controller (SSC)
by [10]. The control inputs which stabilize the position of Rj relative to ~ are [11]
86
y
where
'+x
Figure 7.1: Definitions for the SBC controller (from [10]).
Sij cos 'Yij  lij (bij +Wi) sin 'Yij + Vi cos (()i  (}j)
Sij sin 'Yij + lij (hij +Wi) COS 'Yij + Vi sin (()i  ()j)
dj
(7.1)
'Yij 'I/J .. + (J.  (). tJ t J
ij k1(itj  l