Coho: A Verification Tool for Circuit Verification by Reachability Analysis by C h a o Y a n B . S . , Pek ing Univers i ty , 2003 i A T H E S I S S U B M I T T E D I N P A R T I A L F U L F I L L M E N T O F T H E R E Q U I R E M E N T S F O R T H E D E G R E E O F Master of Science i n T H E F A C U L T Y O F G R A D U A T E S T U D I E S (Computer Science) The University of British Columbia Augus t 2006 Ā© C h a o Y a n , 2006 Abstract C O H O is a ver i f icat ion too l for systems modeled by nonl inear o rd inary differen-t ia l equat ions ( O D E s ) . Ver i f icat ion is performed using reachabi l i ty analysis. T h e reachable space is represented by projectagons wh ich are the po lyhedron descr ibed by their pro ject ion onto two d imensional subspace. C O H O approx imates nonl inear O D E s w i t h l inear dif ferential inclusions to compute the reachable state. We re-engineer the C O H O system to provide a robust imp lementa t ion w i t h a clear in -terface and wel l -organized st ructure. W e demonstrate the soundness and robustness of our approach by app ly ing it to several examples, inc lud ing a three d imensional , non- l inear systems for wh ich previous version of C O H O fai led due to numer ica l sta-b i l i ty problems. T h e correctness of C O H O s t rongly depends on the accuracy, eff iciency and robust-ness of the l inear p rogram solver that is used throughout the analysis. O u r C O H O l inear p rogram solver is implemented by integrat ing the S imp lex a lgor i thm w i th interval a r i thmet ic based on the framework set up by L a z a [LazOl] . For h igh ly i l l -condi t ioned problems, an approx imat ion that is very close to the op t ima l result is computed by our a lgor i thm. T h i s results i n a very sma l l error bound that allows C O H O to successful ly ver i fy interest ing examples. T h i s thesis also presents an a lgor i thm to solve the prob lem of pro ject ing the feasible region of a l inear program onto two d imensional subspaces. T h i s a lgor i thm uses the C O H O l inear p rogram solver for efficiency and accuracy. W e der ive an analy t ica l upper bound for the error and present exper imenta l results to show that the errors are negl igible i n pract ice. i i Contents Abstract ii Contents iii List of Tables vi List of Figures vii Acknowledgments ix Dedication x 1 Introduction 1 1.1 Mo t i va t i on 1 1.2 Con t r i bu t i on 3 1.3 Thes is Organ iza t ion 4 2 Background and Related Work 7 2.1 F o r m a l Ver i f icat ion 7 2.2 Ver i f icat ion as Reachab i l i t y 10 2.3 Pro jec tagons 12 2.4 C O H O 14 3 Architecture of the Coho System 23 i i i 3.1 Java Componen t 24 3.2 N u m b e r Package : 25 3.3 M a t r i x Package 29 4 Coho Linear Program Solver 32 4.1 L inea r P r o g r a m m i n g and the S implex A l g o r i t h m 33 4.1.1 C O H O L inear P rog ramming 33 4.1.2 P r i m a l - D u a l M e t h o d 34 4.1.3 S imp lex A l g o r i t h m 35 4.2 C o m b i n a t i o n of S implex and Interval C o m p u t a t i o n 36 4.2.1 Amb iguous P i vo t i ng and C y c l i n g 38 4.2.2 Poss ib le O p t i m a l Resu l t 40 4.2.3 In i t ia l Feasible Bas is 40 4.3 Eff ic ient L inea r System Solver 42 4.4 I l l -Cond i t ioned Basis and Excep t ion H a n d l i n g 44 4.4.1 Cond i t i on ing 44 4.4.2 Use of Non-op t ima l Bases 47 4.5 Implementat ion 50 5 Projection 53 5.1 A l g o r i t h m 53 5.2 Implementat ion 57 5.3 Conver t ing f rom End-o f -S tep to Beginning-of-Step Pro jec t i on . . . . 61 5.4 E r r o r Ana l ys i s 64 6 Experiments 69 6.1 A p p l i c a t i o n of C O H O Ver i f icat ion System 69 6.2 Expe r imen ta l Resu l t for L P solver and L P project 74 iv 7 Conclusion 80 7.1 W h a t has been Accomp l i shed 80 7.2 Suggestions for Fur ther Research 81 Bibliography 86 v List of Tables 2.1 Compar i son of H y b r i d System Ver i f icat ion Tools 9 3.1 Func t ions of CohoDoub le ln te rva l 29 v i List of Figures 2.1 A Three D imens iona l "Pro jec tagon" 13 2.2 A T i m e Step of C O H O 16 2.3 T h e Crea t ion of a B loa ted L inear P r o g r a m f rom Pro jec tagon . . . . 17 3.1 T h e Arch i tec tu re of C O H O System 23 3.2 Package Hiearchy of C O H O Sys tem 25 3.3 C lass H iearchy of N u m b e r Package 26 3.4 C lass H iearchy of M a t r i x Package 31 4.1 P i v o t Select ion of C O H O L inear P r o g r a m Solver 39 4.2 C O H O L inea r System Solver 44 4.3 T y p e s of O p t i m a l 2D Vert ices 49 4.4 H i g h l y Ob tuse O p t i m a l 3D Vert ices 50 4.5 C lass H iearchy of C O H O L inear P r o g r a m Solver 51 5.1 P ro jec t i on A l g o r i t h m Us ing L inear P r o g r a m m i n g 54 5.2 F i n d the N o r m a l of Cur ren t Edge Us ing O p t i m a l Bas is of Prev ious L inear P r o g r a m 55 5.3 Overes t imat ion of Pro jec t ion Po lygon W h e n an E d g e is Sk ipped . . 57 5.4 T h e Representat ion of Op t im iza t i on D i rec t ion 57 5.5 F i n d the In i t ia l Op t im i za t i on D i rec t ion for P ro jec t i on A l g o r i t h m . . 58 5.6 Back t rack ing to Reduce E r ro r W h e n J u m p to Nex t Quadran t . . . . 60 v i i 5.7 Spec ia l Case of J u m p i n g T w o Quadrants Cont inuous ly 61 5.8 Convers ion between End-o f -S tep Pro jec t ion and Beginning-of -Step P ro jec t i on 63 5.9 P ro jec t i on E r r o r In t roduced by C O H O L inear P r o g r a m Solver . . . . 65 5.10 Pro jec t i on E r r o r In t roduced by Ignor ing Sma l l Ang le 66 5.11 P ro jec t i on E r r o r In t roduced by Sk ipp ing Short Edge 66 5.12 Pro jec t i on E r r o r In t roduced by J u m p i n g to Nex t Quadran t 67 6.1 T h e Resu l t of a T w o D imens iona l L inear M o d e l E x a m p l e : S ink . . . 70 6.2 T h e Resu l t of 3 V d P E x a m p l e : x-y plane 73 6.3 T h e Resu l t of 3 V d P Examp le : y-z plane 73 6.4 E r r o r D i s t r i bu t i on of O p t i m a l Resul t of C O H O L inea r P r o g r a m Solver 75 6.5 D i s t r i bu t i on N u m b e r of Branches of a P i vo t 76 6.6 D i s t r i bu t i on of Cond i t i on N u m b e r for I l l -Cond i t ioned P rob lems . . . 77 6.7 N u m b e r of L inear Programs and Except ions per Step 78 v i i i Acknowledgments W i t h o u t encouragement f rom M a r k Greenstreent, i t wou ld be impossib le to wr i te th is thesis. It is dif f icult to express my grat i tude for his extensive suppor t , d iscussion and endless intel l igent ideas. I a m grateful to Ian M i t che l l , C h e n Grei f , and M i c h a e l Fr ied lander for their cont r ibu t ion of t ime and ideas to my work. I wou ld also l ike to thank M a r i u s L a z a , Ian M i t c h e l l and many others I do not know their names for their previous hard work on C O H O . W i t h o u t their work, it is impossib le to complete C O H O by myself. Spec ia l thanks to m y parents and yu l iang for too many reasons to say. I wou ld l ike to thank my fr iends, Suwen Y a n g , X i u s h a n Feng , Su l ing Y a n g , X u n S u n , W e i L i , X i a n g X u a n , Fei M a , Shu 'an W a n g , J ie Zhao, J u l i a C h e n , C h e n Y a n g and many others, for enjoying the happy life together. C H A O Y A N The University of British Columbia August 2006 To my love, Yu l i ang . x Chapter 1 Introduction 1.1 Motivation W i t h the increasing of the complex i ty of c i rcui t design, i t is more and more impor -tant to veri fy that a c i rcui t satisfies i ts high-level speci f icat ion automat ica l ly and formal ly. P ro to t yp ing and s imulat ion are two means of va l ida t ion . However, neither gives fu l l assurance that a l l possible si tuat ions have been checked and the system always works as expected. Fur thermore, compared w i t h increasing of design com-plexity, the percentage of prob lem cases covered by manua l l y generated test suites or pseudo- random funct ional test ing is decl in ing. Therefore, fo rmal ver i f icat ion, a mathemat ica l approach that exhaust ively proves funct ions propert ies, has at t racted more and more at tent ion. For d ig i ta l c i rcui ts or hybr id systems, wh ich conta in bo th discrete compo-nents and cont inuous components, the veri f icat ion prob lem can be formulated as reachabi l i ty analysis. T h e system is modeled by some mathemat ica l mode l , such as ord inary dif ferential equat ions ( O D E s ) . Thus , ver i f icat ion amounts to comput ing the set of reachable states and checking the propert ies on such states. C O H O , a ver i f icat ion too l for c ircui t design, performs reachabi l i ty analysis. It models c i rcui ts by O D E s and represents reachable regions by their project ion onto two-d imensional subspaces. General ly, the nonl inear O D E mode l does not 1 have closed fo rm solut ions, and approx imat ion techniques must be used. C O H O computes an over approx imat ion of the reachable space such that it might fai l to veri fy a correct system, but it never mistakenly verifies an incorrect design. To compute the evolut ion of reachable states, C O H O solves a large number of l inear programs. T h e accuracy and efficiency of the l inear p rogram solver plays a great role for the correctness and performance of C O H O . L inea r p rogramming is a wel l -know prob lem w i t h classical solut ions such as the S imp lex a lgor i thm. E x -p lo i t ing the specia l s t ructure of l inear programs that arise f rom C O H O , L a z a [LazOl] presented an efficient 0(n) l inear system solver for S imp lex a lgor i thm such that the op t ima l result is accurate and the l inear program solver is efficient. F l oa t i ng point computa t ion does not guarantee that the op t ima l value for the l inear p rogram is a lways over approx imated as requi red. A n under approx imat ion of the error might make it possible that C O H O verifies an incorrect design. Fur ther-more, the soundness of C O H O replies heavi ly on the accuracy of the op t ima l result of l inear programs. Therefore the l inear program solver is required to produce solu-t ions as accurate ly as possible, especial ly for the bad ly condi t ioned problems that often arise i n the C O H O ' s analysis. Otherwise, the large over -approx imat ion error prevents C O H O f rom ver i fy ing correct systems. These problems are solved by implement ing the S imp lex a lgor i thm w i t h in -terval a r i thmet ic wh ich provides bo th the result and its error b o u n d . Interval com-puta t ion guarantees over approx imat ion of the f inal result . I l l -condi t ioned problems are easi ly detected by the interval representat ion. B u i l d i n g on the interval ar i th-met ic methods, we present a new a lgor i thm for f ind ing a good approx imat ion of the op t ima l result . Us ing the C O H O l inear program solver, we also developed an algo-r i t h m to compute the project ion of a po lyhedron onto a two-d imensional subspace. We demonstrate these capabi l i t ies of these new algor i thms by app ly ing them to the analysis of a s imple, non- l inear system wh ich could not be analysed by pr ior versions of C O H O . 2 ā¢ 1.2 Contribution T h i s research focused on improv ing and implement ing a robust and efficient l inear program solver based on previous work [LazOl] and developing a new project ion a lgor i thm for C O H O . T h i s improvement and development are complemented by a re-engineering of the C O H O source code to produce a robust imp lementa t ion of the C O H O system. T h e ma in cont r ibut ions of th is thesis include: 1. T h e re-engineering of the C O H O software. ā¢ W e defined general purpose interfaces for numer ic types and matr ices. These prov ided an easi ly extended framework for in tegrat ing interval a r i thmet ic into C O H O and implement ing the l inear program solver and pro ject ion a lgor i thms. ā¢ Interval ar i thmet ic is integrated into bo th l inear program solver and pro-ject ion a lgor i thms. ā¢ Re-engineer ing of the l inear program solver to improve its modular i ty . T h e var ious phases of the S implex a lgor i thm are c lear ly separated, and we make systemat ic use of exceptions for hand l ing very bad ly condi t ioned bases. 2. Improvement and re implementat ion of the l inear p rogram solver us ing interval computa t ion ā¢ T h e solver provides bo th upper bounds and lower bounds wh ich guarantee overapprox imat ion as required by C O H O , and the error bound is quite smal l . ā¢ W e developed a much more efficient a lgor i thm to est imate the condi t ion number of a basis than the previous method. 3 ā¢ We implemented the method proposed in [LazOl] for d iscard ing (nearly) redundant constraints to handle very bad ly condi t ioned bases. ā¢ We took special care to el iminate a l l tests for f loat ing-point equal i ty f rom these algor i thms. T h i s corrected several errors that remained in the pre-v ious implementat ion. 3. A n implementa t ion of a special ized a lgor i thm for pro ject ing the feasible poly-hedron onto a two-dimensional subspace. ā¢ U s i n g our C O H O solver, the project ion method is efficient. ā¢ T h e a lgor i thm over approximates the pro ject ion po lygon w i t h smal l error. 4. A p p l y the revised C O H O system to two examples: ā¢ A s imple, s ink w i th l inear dynamics f rom [DM98] . W e show that our implementa t ion of C O H O provides much t ighter bounds than those pre-v ious ly reported for other tools. ā¢ A van der P o l osci l lator w i th three state var iables. T h i s uses the pro-ject ion capabi l i ty of C O H O . Fur thermore, the previous version of C O H O was unable to ver i fy th is example due to numer ica l s tab i l i ty issues i n the l inear program solver that we have now overcome. 1.3 Thesis Organization T h i s thesis consists of seven chapters as descr ibed below: ā¢ Chap te r 2 int roduces the C O H O system in wh ich the l inear p rogram solver and pro ject ion a lgor i thm are appl ied. F i r s t , reachabl i ty analysis and its appl icat ion to ver i f icat ion, especial ly c i rcui t ver i f icat ion are in t roduced. T h e n the concept of represent ing reachable space by project ions is descr ibed briefly. F ina l l y , the a lgor i thm and structure of the C O H O ver i f icat ion system is presented. 4 ā¢ Chap te r 3 presents the architecture of the C O H O system. T h e par t i t ion a M a t l a b subsystem and a Java subsystem is exp la ined f irst, fol lowed by the descr ip t ion of packages that we wrote for in terval a r i thmet ic and matr ices. ā¢ Chap te r 4 presents the implementat ion of the C O H O l inear program solver. We first give a br ief in t roduct ion to l inear p rogramming and the specia l s t ructure of the l inear programs that arise i n C O H O . T h e n the problems of in tegrat ing S imp lex w i t h interval ar i thmet ic are presented, fol lowed w i t h our solut ions. To improve the speed of C O H O solver, Laza ' s l inear system solver is employed and improved. T h e a lgor i thm to f ind a good approx imat ion of the op t ima l cost when the op t ima l basis is bad ly condi t ioned is also presented w i t h geometric explanat ions. T h e chapter ends w i th a descr ipt ion of some implementat ion issues. ā¢ Chap te r 5 presents our a lgor i thm for the pro ject ion p rob lem. T h e basic idea is that we can use the feasibi l i ty of the current basis to f ind the next edge. We then describe implementat ion problems and solut ions. M o v i n g the reachable region forward in t ime according to its differential inc lus ion mode l produces a l inear p rogram w i t h a different structure than we have assumed so far. For tu -nately, th is end-of-t ime-step l inear program can be expressed as the product of a l inear p rogram w i t h the special C O H O s t ructure and an easi ly invert ib le m a -t r ix . W e then extend our project ion a lgor i thm to work w i t h l inear programs of th is form. F ina l l y , the error of our project ion a lgor i thm is ana lyzed, and an upper bound on the error is obta ined. ā¢ Chap te r 6 presents exper imenta l result of C O H O system. W e first describe two examples and present the results produced by the C O H O ver i f icat ion tool . T h e n the exper imenta l results for C O H O l inear p rogram solver and project ion a lgor i thm are examined wh ich demonstrate the robustness of our a lgor i thms. 5 ā¢ A chapter of conclusions and research topics i n the future completes th is thesis. 6 Chapter 2 Background and Related Work T h i s chapter describes related research in hybr id system ver i f icat ion and pr ior work on the C O H O system. Sect ion 2 . 1 summarizes exist ing approaches for the formal ver i f icat ion of hyb r i d systems. Sections 2 . 2 th rough 2 . 4 descr ibe the C O H O system. Sect ion 2 . 2 describes how C O H O formulates the ver i f icat ion prob lem for c i rcui ts and hyb r id systems as a reachabi l i ty prob lem for systems modeled by differential inc lusions. A cr i t ica l issue i n such analysis is finding a t ractable representat ion for h igh d imens iona l objects. C O H O uses "pro jectagons" , where the reachable space is represented by i ts pro ject ion onto two d imensional subspaces. Pro jectagons are in t roduced in section 2 . 3 . T h e n , Sect ion 2 . 4 describes how a lgor i thm that C O H O uses to compute the reachable space for a c ircui t or hyb r id system. 2.1 Formal Verification Forma l ver i f icat ion endeavors to prove or disprove that a system satisfies a cer-ta in formal speci f icat ion or proper ty using formal methods f rom mathemat ics . T h e design might be a hardware system, a software system, a network pro toco l , an air-plane, etc. T o veri fy such systems, formal , abstract ma themat i ca l models must be constructed for bo th the system and the speci f icat ion. C o m m o n l y used mathemat-ica l models inc lude f inite state machine ( F S M ) , label led t rans i t ion systems, Pe t r i 7 nets, etc. Recent ly , fo rmal ver i f icat ion for discrete systems have gained indust r ia l suc-cess, especial ly i n the areas of d ig i ta l c i rcui ts and communica t ion protocols. E x t e n d -ing formal ver i f icat ion techniques to systems w i th cont inuous dynamics has been a par t i cu la r ly chal lenging research f ield. Hybrid Systems are dynamica l systems w i t h bo th discrete and cont inuous behaviours. For example, an embedded system, a flight control ler and a c i rcui t w i t h bo th d ig i ta l and analog singals are al l hybr id systems. Ver i fy ing hybr id systems formal ly is increasingly impor tan t because of the growing prevalence of computers as control lers i n safety cr i t ica l systems and the impor tance of analog c i rcui t effects i n deep-submicron integrated circui ts. A common ly used formal mode l for hybr id systems is the hyb r id automa-ton [Hen96]. Ver i f icat ion problems for most classes of systems w i t h nonl inear con-t inuous dynamics are undecidable. A common approach to ver i f icat ion of such systems thus involves some type of approximate reachabi l i ty ca lcu la t ion. M a n y ver i f icat ion tools have been developed for hyb r id systems. Table 2.1 summar izes four of the most wide ly used tools that are closely related to C O H O . HyTech [ H H W T 9 7 ] [ H P W T 0 1 ] is a symbol ic mode l checker for linear hybrid au-tomata [Hen96]. T h e hybr id system is modeled as a col lect ion of au tomata w i t h discrete and cont inuous components, and the propert ies to ver i fy are ex-pressed us ing tempora l logic formulas. H y T e c h can check safety propert ies; if the ver i f icat ion fai ls, H y T e c h generates a diagnost ic error trace. Ano the r key feature of H y T e c h is i ts abi l i ty for parametr ic analysis to determine necessary and sufficient constraints on the parameters under wh ich the system is safe. d/dt [ADM02] is a too l for reachabi l i ty analysis of cont inuous or hyb r id systems us ing l inear differential inclusions of the form of dx/dt = Ax + Bu, where u is an external inpu t tak ing values in a bounded convex po lyhedron, d/dt represents the reachable sets as non-convex or thogonal po lyhedra [BMP99 ] , 8 Table 2.1: Compar i son of H y b r i d System Ver i f icat ion Too ls Too l H y T e c h d / d t C h e c k M a t e P H A V e r Model Linear hybrid au-tomata Hybrid automata with linear differ-ential inclusion Polyhedral in-variant hybrid automata Linear hybrid "au-tomata C T L specifi-cation yes no A C T L no Representation of Reachable Set Hyper rectangular Orthogonal poly-hedron Polyhedron defined by Cx < d Polyhedron de-fined by a set of constraints Approximation Technique 1) clock transla-t ion 2) rate trans-lation 3) linear phase-protrait ap-proximation face lift-ing [DM98] 1) quotient tran-sition systems (QTS) 2) flowpipe approximation On-the-Fly over-approximation of piecewise affine dynamics Counter E x -ample yes no yes no i.e. f ini te unions of fu l l -d imensional hyper-rectangles, and approximates the reachable state using numer ica l integrat ion and po lyhedra l approx imat ion , d / d t can be used to compute reachable sets, ver i fy safety propert ies of systems, and synthesize safety propert ies of swi tch ing control lers based on the maximal invariant set. CheckMate [SR + 00 ] is a M a t l a b based too l for mode l ing , s imu la t ing and ver i fy ing propert ies of a class of hybr id systems: threshold-event-driven hybrid systems. T h e system is represented as a polyhedral invariant hybrid automaton (PIHA), wi th three types of dynamics : clock dynamics (dx/dt = c, l inear dynamics (dx/dt = Ax + b) and nonl inear dynamics (dx/dt ā f(x)). Therefore, ver i-f icat ion can be performed d i rect ly on a mode l of the hyb r i d system wi thout const ruct ing a l inear approx imat ion [SKOO], as H y T e c h and d / d t do. The system speci f icat ion is expressed as an A C T L formula, a restr ic t ion of C o m -puta t iona l Tree Log ic ( C T L ) [Eme97]. C h e c k M a t e can per fo rm two k inds of ver i f icat ions: a quick ver i f icat ion and a complete ver i f icat ion. W i t h the quick 9 ver i f icat ion, on ly the vertices of the po lyhedron defined by the in i t ia l condi-t ions are veri f ied. T h e advantages of the quick ver i f icat ion is that it is fast and can provide a rough idea about the f inal result of the complete ver i f icat ion. If the quick ver i f icat ion finds an error, then the design violates the specif icat ion and should be corrected. However, the quick ver i f icat ion may fai l to find some errors. T h u s i f the quick veri f icat ion succeeds, then a complete ver i f icat ion can be per formed to make sure that the design is correct. PHAVer [Fre05] is a ver i f icat ion tool for l inear hybr id au tomata that supports compos i t iona l and assume-guarantee reasoning. It uses the P a r m a Po l yhed ra L i b r a r y that suppor ts arb i t rary large number, thus it provides exact, robust ar i thmet ic . T h e system mode l is based on affine dynamics that are handled by on-the-f ly overapprox imat ion. Table 3.1 compares some impor tant features of these tools. There is an-other hyb r id system ver i f icat ion tools that deserves some at tent ion: VeriShift. In contrast to al l tools descr ibed above, it approximates the reachable region using el-l ipsoids [BTOO]. T h e advantages of el l ipsoidal approx imat ion are a reduced memory requirement and a po lynomia l complex i ty of e l l ipsoidal operat ions. 2.2 Verification as Reachability T o ver i fy a sys tem, the fo rmal mode l is checked for correctness w i t h respect to the formal ly expressed speci f icat ion. T h e entire reachable state space of the sys-tem shou ld be explored. Represent ing and comput ing these reachable states is the essential step for many ver i f icat ion tasks. G i v e n a ?i-dimensional region A G Rn, the reachabl i ty p rob lem is to compute the (forward) reachable state space that contains al l trajectories that start in A according the mode l of the system, either dur ing some t ime interval [0,Ā£ eāa] o r for al l t ime. Genera l ly , the mode l can be represented by an ord inary differential 10 inc lus ion: xeF(x,t) (2.1) where x ā¬ Rn. T h e inc lus ion i n the mode l allows uncertainty. M a n y ver i f icat ion problems can be solved as reachabi l i ty analysis problems. For example, the safety proper ty ver i f icat ion prob lem is to explore the entire reach-able states space of system and val idate the proper ty on a l l the reachable states. Reachab i l i t y analysis can be used to veri fy whether a c i rcu i t correct ly imple-ments i ts speci f icat ion. General ly , a c i rcui t can be modeled by ord inary differential inclusions (ODIs ) , and the in i t ia l state and constraints determine the in i t i a l region i n the state space. T h e n , reachable sets are computed by advancing the in i t ia l region accord ing to the O D I mode l . A f ter explor ing a l l the reachable states, the c i rcui t ver-i f icat ion prob lem becomes a prob lem of checking whether the speci f icat ion is va l id for a l l reachable states. Unfor tunate ly , closed form solut ions do not exist for general O D I models, except for some special cases such as l inear models. T h u s , an approx imat ion based approach must be appl ied for reachabi l i ty analysis. For ver i f icat ion problems, the reachable state space should be over approx imated. It is sound that incorrect sys-tems are never falsely veri f ied, a l though the over approx imat ion might make it fai l to ver i fy a correct system. In C O H O the reachable space at each t ime step is par-t ioned into smal l pieces, and O D I model for each piece is approx imated by a l inear differential inc lus ion: F(x, t) C A ā¢ x + b + cube(u) (2.2) where cube(u) is a hypercube, i.e., the Car tes ian product of n intervals, wh ich represents an error bound for the approx imat ion. 11 2.3 Projectagons T h e representat ion of the reachable space strongly affects the efficiently and accuracy of a reachabi l i ty analysis a lgor i thm. For example, approx imat ing the reachable space by the m i n i m u m hyper-rectangle that contains it is qui te s imple. T h e number of vert ices/ faces of hyper-rectangles increases l inear ly w i t h the d imension, and operat ions on hyperrectangles are also qui te efficient. However, hyperrectangles produce an unacceptab ly large over approx imat ion of the reachable space for most appl icat ions. A t the other ex-treme, general nonconvex h igh-d imensional po lyhedra can represent arbi t rary, h igh-d imensional objects w i t h relat ively smal l error. However, man ipu la t ions of general n d imens ional po lyhedra typ ica l ly have t ime and space complexi t ies w i t h exponents of n or n / 2 . In C O H O , reachable sets are represented as projectagons. For the purposes of th is thesis, a projectagon is the h igh d imensional (and potent ia l ly nonconvex) bounded po lygon formed by the intersect ion of a col lect ion of pr isms. E a c h p r ism is unbounded i n a l l but two dimensions, and i n those two dimensions the cross-section of the p r i sm is a bounded po lygon. B y choosing th is representat ion we need on ly t rack the two-dimensional cross-sections of the pr isms rather than any ful l d imens ional object. We take advantage of this s impl i fy ing proper ty i n bo th the computa t iona l geometry and numer ica l a lgor i thms. ' For example, F igu re 2.1 shows how a three-dimensional object (the "anv i l " ) can be represented by i ts pro ject ion ā¢ onto the xy, yz, and xz planes. T h e pro ject ion polygons are not required to be convex; thus, non-convex, h igh-d imensional objects can be represented by projectagon. T h e h igh-d imensional object represented by a projectagon is the largest set of points that satisfies the constraints of each pro ject ion. It can be obta ined f rom its project ions by back-pro ject ing each pro ject ion po lygon into a pr ism in Rn and comput ing the intersect ion of these pr isms. T h i s can lead to an overapprox imat ion of the po lyhedron. Thus , 12 not a l l po lyhedra can be represented exact ly using projectagons. For example, the projectagon representat ion of an object w i th one or more indentat ions on its faces is an pro jectagon that has these indentat ions fi l led i n . Pro jectagons offer several advantages. Because projectagons are based on two-d imensional polygons, they can be efficiently man ipu la ted using wel l -known al-gor i thms f rom computa t iona l geometry [PS85]. Ignor ing degeneracies, faces of the object represented by a projectagon correspond to edges of i ts pro ject ion polygons. A s s u m i n g that the mode l is l inear and has finite derivat ives, the ex t remal trajecto-ries are those emanat ing f rom faces. The projectagon representat ion allows many operat ions on faces inc lud ing calculat ions of their t rajector ies to be carr ied out as s imple operat ions on po lygon edges. M a n y other representations have been used in reachabl i t i ty analysis tools. Bournez and M a l e r [BMP99] investigates orthogonal polyhedra, wh ich par t i t ion poly-hedra into f inite unions of fu l l d imensional hyper-rectangles. T h i s representat ion is 13 canonical for a l l po lyhedra , convex or non-convex, i n any d imension. Geomet r ic op-erat ions on or thogonal computat ions are efficient. T h e y focus on a specia l case that they cal l griddy polyhedra where the vertices of the po lyhedron are constra ined to lie on evenly spaced gr idpoints i n Rd. [SK03] presents a method that approximates reachable sets us ing oriented rectangular hulls ( O R H s ) . A n O R H is a hyperrectan-gle whose or ientat ion is determined by the reachable set rather t han being forced to be para l led to the coordinate axes of the model . A preferred or ientat ion is derived f rom the s ingular value decomposi t ion ( S V D ) of sample covariance matr ices for sets of reachable states. T h e orientat ions keep the over-approx imat ion of the reachable sets sma l l in most cases w i t h a complex i ty of low po lynomia l order w i t h respect to the d imens ion of the cont inuous reachable state space. A s ment ioned in section 2.1, e l l ipsol idal techniques are also at tract ive because they use on ly 0 ( n 2 ) space for n dimensional object , and the complex i ty of e l l ipsoidal operat ions is also po lynomia l rather t han exponent ia l for general po lyhedra l operat ions. However, the d isadvan-tage is that the intersect ion or un ion of el l ipsoids are no longer el l ipsoids, and thus approx imat ions must be performed. 2.4 Coho C O H O [GM98 , G M 9 9 ] is a reachabl i ty analysis too l for systems modeled by non-l inear o rd inary dif ferential equat ions. C O H O uses projectagons to represent reach-able regions, and over approximates the ord inary dif ferential equat ions mode l by l inear dif ferential inc lus ion, as descr ibed i n section 2.2 and 2.3. A s descr ibed i n the previous sect ion, the projectagon representat ion allows many operat ions on faces inc lud ing calculat ions of their t rajector ies to be carr ied out as s imple operat ions on po lygon edges. Therefore, C O H O advances the pro-jectagon of a reachable set by consider ing each face of the pro jectagon in tu rn . C O H O moves each edge/face forward in t ime w i t h an over approx imat ion obta ined by a s imple E u l e r integrator [RLB01] . A l t h o u g h an E u l e r integrator is less accurate 14 t han higher-order methods, its s impl ic i ty allows us to guarantee that our results are overapproximat ions as desired. A s i l lus t ra ted in F igure 2.2 f rom [ L a z O l ] , C O H O advances the reachable region by per forming the fol lowing operat ions at each t ime step (see [GM99] for detai ls) : 1. Bloat the projectagon Le t P be a projectagon. To compute al l points reachable f rom P i n the next t imes step, we first compute H, an possib ly over approx imat ion of the convex hu l l of P, and then B, a bloated version of H. H is s imp ly the projectagon corresponding whose project ion polygons are the convex hul ls of the projec-t ion polygons of P; i n other words, by approx imate the convex hu l l of P by comput ing the hul ls of the each of its pro ject ion polygons independent ly. We can describe H as a l inear program: Ax < c. W e normal ize A such that each row of A has a I2 no rm of 1. Now, the l inear p rogram for B is Ax < c + dbioat, where duoat is the b loat ing distance. F igure 2.3 i l lustrates this process. In the fol lowing, we refer to a project ion polygon of a projectagon as a slice. 2. Linearize the model for each edge For each edge, e, C O H O creates a rectangle, Re whose edges are either paral le l or perpend icu lar to e such that Re contains e w i t h a separat ion of duoat. Clear ly , Re can be represented by a l inear p rogram; thus, we w i l l wr i te Re to denote th is l inear program in the fol lowing. Re is used to compute the step size such that trajectories wh ich start on the face corresponding to e w i l l rema in in the' intersect ion of Re and B du r ing the current t ime step. C O H O now uses l inear ized, over approx imat ion of the non- l inear mode l : .F(x,t) C A - x + b + cube(u) (2.3) where cube(u) = {z \ V i . < v*}. T h e cube(u) te rm allows the nonl inear re lat ion to be approx imated by a l inear one plus an error te rm. T h i s error 15 advance F(x) C Ax + b + eube(u) project / " -1 .2 0 -2.1 _ X 2.0 ' ^advance = 1 0 -10.3 < 10.8 0 - 5 1.4 y -3~7~ I 0 1 3.2 z oi.9 intersect assemble projections] 7% V \ \ \ \ \ \ bloat each face advance dach slice \ / ā¢-' P bloat union arid simplify .clip slices bnext x e F(x) F(x) C Ax + b + cube(u), V x ā¬ euoat Figure 2.2: A T i m e Step of C O H O 16 F igu re 2.3: T h e Crea t ion of a B loa ted L inear P r o g r a m f rom Pro jec tagon [LazOl , F igure 2.2] 17 t e rm is usua l ly very smal l because the variables are restr ic ted to the smal l rectangle R; thus, the nonl inear O D E can be approx imated by l inear one w i th a smal l error bound . If the edge is too long to ob ta in a good approx imat ion , C O H O spl i ts it in to several short edges. T h e l inear izat ion of the mode l is performed by code prov ided by the user. C O H O cal ls a user prov ided M a t l a b funct ion, createmodel, passing it the l inear p rogram for the intersect ion of Re and B. T h e createmodel funct ion can invoke the l inear p rogram solver to obta in bounds on the var iables (and their l inear combinat ions) . Based on these bounds, createmodel computes values for the ma t r i x A, and vectors b and u for equat ion 2 . 3 . These are re turned to C O H O for fur ther use i n comput ing the t ime step. It the user provides a createmodel funct ion that is incorrect , then C O H O can produce inva l id results. 3 . Compute the step size For each face, given the l inear ized mode l and the b loa t ing distance, dbioat, C O H O can determine a step size such that al l t rajector ies s tar t ing f rom this face are guaranteed to stay inside the region B f l Re du r ing the current t ime step. T h e next step size At is the m i n i m u m value of a l l such step size for each face, and the next step t ime tnext equals tcurr + At. Therefore, the reachable space at the end of the t ime step must be contained in B D Re. 4. Advance each slice G i v e n the l inear ized mode l f rom step 2 and the step size f rom step 3 , C O H O computes the t ime-advanced slice w i t h the fol lowing operat ions: (a) Advance each face C O H O advances each slice by advancing each edge of the project ion po lygon for that edge (i.e. each face of the projectagon). Because the or ientat ion of the face may rotate dur ing the t ime step, i ts pro ject ion at the end of the t ime step can be a po lygon, rather t han an edge. Thus , we 18 refer to the post- image of an edge as an edge polygon. W e first exp la in how C O H O computes the edge polygons. We then descr ibe how C O H O combines these edge polygons to obta in a t ime advanced slice. i. Bloat the edge R eca l l that for each face, there is a corresponding edge of a project ion po lygon. Le t e be such an edge. C O H O constructs the rectangle Re as descr ibed in step 2. i i . Intersect the bloated edge with original polygon T h e intersect ion of Re and the or ig ina l po lygon for th is slice is computed. T h e result po lygon contains the current face and a smal l inward region. It might be non-convex, i n wh ich case, i ts convex hul l is computed and used in the fol lowing steps. i i i . Advance the intersection according to the linearized model T h e intersect ion f rom the previous sect ion can be descr ibed by the fol lowing constraints: Mx0 > q (2.4) For s impl ic i ty , we first expla in how the edge is advanced neglect ing the approx imat ion error, cube(u). We then descr ibe how C O H O in -corporates this error term. W i t h this s impl i f ica t ion, the l inear mode l f rom step 2 becomes: x = Ax + b (2.5) T h e advanced region can be computed by integrat ing. Le t xe be the pos i t ion at t ime tnext of a point whose pos i t ion is x$ at t ime tcurr. B y integrat ing equat ion 2.5 for the t ime per iod \tcurr, tnext], we obta in : xe = e A A t x 0 + (eAAt - I)A~lb (2.6) T h e combinat ion of equat ion 2.4 and equat ion 2.6 produces the con-19 straints for the advanced face Fadvance. MExe>qe (2.7) where E = e~AAt and qe = q + M(I - e ' ^ A ' 1 A t. O f course, the actual mode l is l inear dif ferential inc lus ion, error bound must be considered for the advanced face. C O H O does this by t reat ing the error te rm as an input to the sys tem. A t any t ime instant, the worst-case such input w i l l correspond to one of the cor-ners of the error cube. A l t hough the worst-case corner may change over the course of a t ime step, C O H O makes the s imp l i f y ing assump-t ion that t ime steps are smal l , and that the same corner w i l l be the worst-case one throughout the t ime step. W i t h th is assumpt ion, C O H O treats each d imension of the error cube separately to approx i -mates the worst-case error due to non- l inear i ty in the mode l . T h i s is the one place where C O H O could potent ia l ly compute a sl ight under approx imat ion of the reachable space. T h i s is an issue that has been present since the or ig inal design of C O H O [GM99] . Address ing this potent ia l unsoundness is beyond the scope of the present thesis and is a topic for future work. iv . Project back the advanced face onto the slice T h e feasible set for Fadvance is a convex po lyhedron . A s descr ibed i n Chap te r 5, C O H O computes the project ion of th is po lyhedron onto the two-dimensional subspace of th is slice. v. Intersect the projection polygon with the bloated slice W e know that the advanced face should be conta ined w i t h i n the b loated region. Therefore, the intersect ion of the pro ject ion polygon and the b loated polygon is computed to get a t ighter approx imat ion , (b) Construct the advanced slice W h e n the t ime-forward images of al l of the edges of the slice have been 20 computed , C O H O takes the un ion of these "edge po lygons" to determine the bounda ry of the polygon at the end of the t ime step. Typ ica l l y , this increases the number of edges in the po lygon. T h e n , C O H O performs overapprox imat ing s impl i f icat ions to keep the number of edges in the po lygon t ractable. 5. Clip advanced slices and create new projectagon U s i n g the steps descr ibed above, C O H O produces a new po lygon for each slice of the pro jectagon. These polygons define a projectagon that contains the reachable space at the end of the t imestep. However, these polygons might not be feasible for each other. Therefore, C O H O c l ips them (project the pro jectahedron onto each two d imensional subspace) to ensure that they are mu tua l l y feasible. T h e projectagon Pnext of the next t ime step is created f rom these c l ipped polygons. A s seen f rom this short descr ipt ion, C O H O makes extensive use of l inear programs when comput ing reachable sets for systems modeled by non- l inear O D E s . For example, step 4(a) i i , 4(a)iv, 4(a)v and 5 need the funct ion to compute project ion of a projectagon onto a two d imensional subspace. T h e pro ject ion funct ion and step 3 use lots of l inear programs. W h e n C O H O was first implemented [GM99] , Greenstreet and M i t che l l found that the l inear programs that arise i n the reachabi l i ty computa t ion are often h igh ly i l l -condi t ioned. W h i l e they were able to analyse a few s imple systems, this i l l -cond i t ion ing prevented further appl icat ion of the approach. M o r e recent ly [LazOl] , the specia l s t ructure of the l inear programs ar is ing i n C O H O was used to obta in efficient and robust solut ions. However, the error was st i l l too large for bad ly condi-t ioned problems, and the error bounds where not made avai lable to the project ion a lgor i thm. T h u s , the pro ject ion step often fai led when appl ied to an i l l -condi t ioned l inear p rogram. These problems mot ivate us to apply interval a r i thmet ic to the l inear pro-21 gram solver and design new a lgor i thm for i l l -condi t ion system hand l ing and proj t ion computa t ion . 22 Chapter 3 Architecture of the Coho System T h e previous chapter inc luded a descr ipt ion of pr ior work on the C O H O ver i f icat ion too l that is the subject of this thesis. T h i s chapter describes the organizat ion of the C O H O software w i t h an emphasis on the packages that were developed or modi f ied i n the current research. approximate the model construct projectahedron advance faces of projecta-hedron at each step probl o solve Solve linear programs Compute projection Gemometry computation result > for Matlab \ Matlab Component FIFO Java Component Figure 3 . 1 : T h e Arch i tec ture of C O H O Sys tem A t the top- level , C O H O has two ma in components: a component wr i t ten in M a t l a b and a component wr i t ten i n Java, as shown in F igure 3 . 1 . T h e M a t l a b com-ponent provides the user interface. In par t icu lar , the mode l for the system being veri f ied is wr i t ten as a M a t l a b module. T h e user interacts w i t h C O H O th rough a set of M a t l a b funct ions. These funct ions suppor t creat ing a mode l , def in ing projectahe-2 3 d ra and comput ing reachable regions. The Java component provides funct ional i ty that is not readi ly available in M a t l a b . In par t icu lar , i t provides operat ions for man ipu la t i ng the pro ject ion polygons and the robust l inear program solver that is the focus of this thesis. T h e Java and M a t l a b components communicate th rough a pai r of named pipes. There is a C program invoked by M a t l a b process that creates these pipes. Op t ions prov ided on the command l ine of the C program al low the content of the streams to be logged. Fur thermore, the communica t ion between the two programs uses p la in A S C I I hence is human readable. T h i s capabi l i ty faci l i tates developing and debugging the Java and M a t l a b funct ions separately. To use the funct ions f rom the Java component , the M a t l a b component wri tes an expression on the m2j (Mat lab- to-Java) p ipe. T h e Java component includes a simple interpreter that evaluates these expressions and wri tes the results back on the j2m (Java- to -Mat lab) p ipe. M a t l a b and Java programs can send numer ica l values either using scientif ic nota t ion or as a hex encoding of the bi ts . T h e later makes the files harder to read but ensures absolutely no loss of precis ion. 3.1 Java Component F igu re 3.2 shows the top level structure of the Java component . O u r robust l inear program solver is implemented i n the Ip packages and uses classes f rom the number and matrix packages. T h e Java component consists of four subsystems: an interpreter for mat lab-java interface, geometry computa t ion packages, in f rast ructure packages, and our l inear p rograming packages. T h e interpreter employs J L e x [BA] and C U P [SHA] for lex ical analysis and pars ing. It interprets commands, f rom the m2j (Mat lab- to-Java) p ipe and calls the corresponding funct ions f rom other packages. T h e o ld version interpreter is used w i t h on ly minor modi f icat ion and it is not descr ibed further here. T h e geometry package is used for polygon operat ion, after each advance step. 24 F igu re 3.2: Package Hiearchy of C O H O Sys tem It make extensive use of the Ben t l ey -O t tman a lgor i thm [B079] [PS85, 271-277]. T h e number and ma t r i x packages are described in the fo l lowing sections. T h e l inear program solver is presented in Chapte r 4. T h e l inear program pro ject ion is descr ibed in Chap te r 5. 3.2 Number Package To implement the l inear program solver, interval ar i thmet ic is employed to t rack computa t ion error. Because round off error is inevi table for f loat ing point com-puta t ion , the error also accumulates and no error bound is prov ided. For C O H O ver i f icat ion sys tem, under-approx imat ion is unacceptable; on the other hand , inter-va l computa t ion is rel iable because it produces bo th results and error bounds , wh ich al low us to guarantee over approx imat ion. A new interface is required to implement interval number . Interval number can not extend the bu i ld - in da ta type Java.lang.Double (also Integer, Boolean, etc) because it is a f inal class. Fur thermore, Java.lang.Number does not provide the ar i th-met ic operat ions that we need to define s tandard ma t r i x operat ions and implement our l inear p rog ramming a lgor i thm. T h e new interface is cal led CohoNumber rather than Number to avoid name confusion. CohoNumber extends Java.lang.Comparable. Several implementat ions of C o h o N u m b e r are prov ided in the number pack-age. T h e inher i tance st ructure of these classes is show in F igu re 3.3. C o h o N u m b e r 25 CohoNumber BasicNumber CohoBoolean Coholnteger CohoLong CohoDouble CohoDoubleLong CohoType BasicType Ā« i n t e r f a c e Ā» CohoNumberl +add() +sub() +mult() +div() +abs() +negate() +sqrt() +max() +min() +random() h ā c H Ā« i n t e r f a c e Ā» Comparable! CohoBoolean r -v: boolean i Coholnteger f I -v: int d Ā« a b s t r a c t Ā» BasicNumber] +create() +convert() +equals() +compareTo() 5 ~ ~ 1 CohoDouble -v: double i 1 CohoLong \ -vi long CohoDoublelnterva - l o : double +hi: double ~) +E() +isDouble() +intersect() +union() +overlap() +contain() +greaterThan() +lessThan() F igu re 3.3: Class Hiearchy of N u m b e r Package 26 provides funct ions for s tandard ar i thmet ic , inc lud ing add, sub, mult, div, abs, max, min, etc. T h i s allows other classes to bu i l d on C o h o N u m b e r and work w i t h any under ly ing representat ion of the data. T h e BasicNumber class defines an abstract number that implements the Co-hoNumber interface. It al lows a new implementat ion of C o h o N u m b e r to provide a m i n i m a l set of ar i thmet ic and compar ison operat ings. P r o m o t i o n is performed automat ica l ly when a b inary operator is appl ied to different k inds of C o h o N u m -bers for extensibi l i ty. T h e order of p romot ion, denned in the BasicType class, is CohoBoolean ā> Coholnteger ā> CohoLong ā> CohoDouble ā> CohoDoublelnterval, wi th possible lost of precis ion when promot ing a CohoLong to a CohoDouble for longer number n w i t h |n | > 2 5 3 . If a new C o h o N u m b e r is added, we do not need to mod i fy the BasicNumber or other CohoNumbers to provide the appropr ia te promo-t ions. A f te r denn ing the promot ion sequence in the BasicType class, the imp lemen-ta t ion is independent. O u r ma in concern is CohoDoublelnterval, the class for . interval numbers. C o -hoDoub le ln te rva l is a representat ion of number i by a range f rom its lower bound to its upper bound . In the fol lowing, we use blodface l ike x for interval values, and x, x for i ts upper and lower bound . CohoDoublelnterval implements the four basic operat ions as: x + y = [x + y, x + y] x - y = \x-y,x-y] x x y = [mm(xy, xy, xy, xy), m a x ( x y , xy, xy,xy)] l ' / x = [l/x,l/x] iixx>0 Obvious ly , these funct ions satisfy the inclusion property [DMn05] : x Ā» y = {xĀ®y\x Ā£ x ,y e y j , for Ā® S { + , - , x , / } (3.1) However, for general interval ar i thmet ic , the lower (upper) bound can not be represented exact ly i n computer . Therefore, round off error must be considered for 27 in terval computa t ion . O u r implementat ion uses the proper ty that the I E E E floating point s tandard requires round ing to the nearest, representable number . Thus , the double precis ion result is w i th in \ulp (unit of least precision) of the exact value. However, the smallest value to be added is a single ulp; thus, the interval is widened by 2ulp du r i ng each operat ion. Several other operat ions are also implemented as: m a x ( x , y ) m in ( x , y ) x U y -x, āx\ [min(|x|, | x | ) ,max( |x | , \x\)] if xx > 0 [0, max( |x | , \x\)} iixx<0 [ m a x ( i , y), max (S , y)] [m in (x , y ) ,m in (S , y)} [min(x ,y ) ,max(x ,y ) ] ] [max(x,y),min(3:,y)}} if max(x,y) < mm(x,y) except ion otherwise [y/x, y/w] if X > 0 except ion otherwise Compar i son operat ions for interval are different because the result can be ambiguous. If two intervals are overlap, the compareTo funct ion throws a Not-ComparablelntervalException. A d d i t i o n a l compar ison operat ions are also prov ided, inc lud ing greaterThan, lessThan, geq, leq, neq etc. Table 3.1 shows the detai ls. T h e I E E E floating point s tandard [Gol91] defines round ing modes inc lud ing round-up and round-down. T h u s , a hand-coded imp lementa t ion tha t used the fu l l capabi l i ty of the hardware wou ld probab ly be faster and provide somewhat t ighter error bounds . T h i s w i l l be descr ibed in greater deta i l i n Chap te r 7. x n y = 28 Table 3.1: Funct ions of CohoDoub le ln te rva l Func t i on Descr ip t ion E x a m p l e g rea te rThan return true if x is clearly greater than y; oth-erwise, return false . [ 2 , 3 ] is greater than [ 0 , 1 ] but not greater than [ 1 , 2 ] l essThan return true if x is clearly less than y; return false otherwise. [ 2 , 3 ] is less than [ 4 , 5 ] but not less than [ 3 , 4 ] geq greater than or equal ' [ 2 , 3 ] is greater than or equal [ 1 . 2 ] leq less than or equal [ 1 , 2 ] is less than or equal [ 2 , 3 ] neq not the same double value [ 1 , 2 ] is not equal with [ 1 , 2 ] same return true if x and y have the same lower and upper bound [ 1 , 2 ] is the same with [ 1 , 2 ] compareTo return 1 if x is clearly greather than y, return - 1 if x is clearly less than y; return 0 if x and y are the same double value, throw exception otherwise [ 1 , 3 ] compareTo [ 2 , 4 ] throws an exception equal return true if x and y are the same double value, return false otherwise [ 1 , 2 ] is not equal with [ 1 , 2 ] , [ 1 , 1 ] equals [ 1 , 1 ] overlap return true if the intersection of x and y is not empty [ 1 , 2 ] and [ 2 , 3 ] overlaps conta in return true if x contains y [ 2 , 4 ] contains [3,4] 3.3 Matrix Package T h e matrix package provides mat r i x computat ion based on the number package. T h e elements of C o h o M a t r i x are of instances of C o h o N u m b e r . T h i s package includes a Matr ix interface and several k inds of matr ices. T h e package st ructure is shown in F igure 3.4. Bas i c ma t r i x operat ions are prov ided, such as get or assign value; ar i thmet ic operat ions inc lud ing transpose, add , sub, mul t , etc; operat ion for rows (columns) and compar ison operat ions. BasicMatr ix implements most funct ions that are independent of the specif ic element da ta type. Mapper is used to implement elementwise operat ions, such as the add i t ion of two matr ices. Reduce is for operat ions that compute a single result f rom al l of the elements, such as the m a x i m u m element of a mat r i x . U s i n g this capabi l i ty, very few methods of B a s i c M a t r i x are abstract. In pract ice, a concrete subclass of 29 B a s i c M a t r i x on ly need to provide methods that assign a value to an element, get the value of an element, and create an object of the element type. T h i s makes it very easy to create new version of M a t r i x . O f course, the subclass can override some of the default imp lementa t ion of methods for efficiency. For example, CohoMatr ix overrides the create funct ion because we do not need to create the zero elements for a C o h o M a t r i x . 30 3 cm' c CD W !^ O 8 s o cr CD Matrix BasicHatrix BooleanMatrix IntegerMatrix DouleMatrix DoublelntervalMatrixj Mapper MapCol MapRow Reduce ReduceCol ReduceRow SparseMatrix Range MatrixValue MatrixError Ā« i n t e r f a c e Ā» Matrix +v() +assign() +arithemetic ops() +compare ops() +col/row ops() + . . .0 Ā« i n t e r f a c e Ā» SparseMatrix +nonZero() +isZero() +nonZeroNumOfRow() +nonZeroNumOTCol() DoublelntervalMatrix +data: CohoDoublelnterval[][] +compareTo() Ā« a b s t r a c t Ā» BasicMatrix -maxReduce: BasicReduce -minReduce: BasicReduce -prodReduce BasicReduce -sumReduce: BasicReduce -binaryMap: Mapper +create() +map() +reduce() + f i l l ( ) 0 CohoMatrix in coho.lp.solver -pos: int [ ] [ ] - isTrans: boolean +getSolution() I BooleanMatrix data: boolean!][] +find() +all() +any() +and() +or() Ā« i n t e r f a c e Ā» Reduce +first() +middle() +last() 5 Ā« i n t e r f a c e Ā» ReduceCol <H +create() Ā« i n t e r f a c e Ā» ReduceRow +create() BasicReduce +f irst : CohoNumber +create() Ā« i n t e r f a c e Ā» Mapper +create() +f() <ā¢' Ā« i n t e r f a c e Ā» MapRow IntegerMatrix ā¢data: int[][ ] +toDouble() Ā« i n t e r f a c e Ā» MapCol DoubleMatrix -data: doublet] [] +toInterval() Chapter 4 Coho Linear Program Solver A s descr ibed i n Chap te r 2, l inear p rogramming plays an impor tan t role i n C O H O . A l l l inear programs are of a special form that we refer to as a C O H O Linear Program. We implemented an efficient and rel iable solver based on the S imp lex a lgor i thm [Dan63] for it. L a z a has set up the framework for C O H O solver. W h a t has been accompl ished by L a z a inc ludes: 1. A l g o r i t h m and implement ion of an O(n) l inear system solver for C O H O ma t r i x [LazOl , pp.32-36]. 2. A l g o r i t h m to deal w i t h uncerta inty and avoidance of cyc l ing [LazOl , pp.87-90]. 3. E r r o r analysis of C O H O solver [LazOl , pp39-81]. W e extend Laza ' s a lgor i thm and re implemented the solver using interval ar i thmet ic . T h e cont r ibu t ion includes: 1. Redes ign the architecture of C O H O solver. 2. Implement the solver based on interval numbers, us ing the Number and Matrix package descr ibed in Chap te r 3.2,3.3. New solver interface and pivot rule are proposed. 32 3. E x t e n d s the l inear system solver for bo th C O H O ma t r i x and C O H O dua l mat r i x , implements a hyb r id solver to reduce error. 4. Dea l i ng w i t h h igh ly i l l -condi t ioned l inear programs; a new method to est imate the cond i t ion number is presented. T h i s chapter is s t ructured as follows: L inear p rog ramming and S imp lex algo-r i t h m is in t roduced in the first section. Sect ion 4.2 describes solut ions for problems caused by interval ar i thmet ic . Laza 's l inear system solver is presented in section 4.3. C O H O l inear programs can be bad ly condi t ioned, so lut ion is presented in section 4.4. A t the last sect ion, implementa t ion issues are descr ibed. 4.1 Linear Programming and the Simplex Algorithm 4.1.1 Coho Linear Programming A L inear P r o g r a m m i n g ( L P ) prob lem is an op t im iza t ion prob lem in wh ich the ob-ject ive funct ion and constraints are al l l inear. A l inear p rogram of s tandard form is m i n c T a ; s.t. Ax = b (4-1) x > 0 T h e co lumn ma t r i x c is cal led the cost vector or optimization direction, whi le the value xopt that minizes cTx is the optimal vertex or optimal point, and cTxopt is the optimal value or optimal cost for the l inear program. T h e constraints define the feasible region of the l inear program. Suppose there are n var iables and m equal i ty constraints, then the feasible region for s tandard l inear p rogram is the non-negative por t ion of a po lyhedron w i th m faces in n-d imensional space. A set of m l inear ly independent columns of ma t r i x A is cal led basis, denoted as B. T h e co lumns in the basis are cal led basic columns whereas others are non-basic. Bas ic co lumns conresponds to basic variables w i t h value of to = A g 1 6 . T h e 33 basic solution for a basis B is defined as 0 itj$B (4.2) to,k if 3 = Bk T h e l inear programs in C O H O express a con junct ion of the half-spaces cor-responding to the project ion polygon edges. These l inear programs can be wr i t ten m in cTx s.t. (4.3) P - E ^ x t b where P-x > b are the constraints such that a; is a point of the projectagon defined in Chap te r 2; c is the cost vector; and E is the forward Eu le r operator for the l inear ized model . Rows of P correspond to project ion po lygon edges; therefore, each row has either one or two non-zero elements. Because the Eu le r steps are relat ively smal l , E is wel l -condi t ioned, and E~l is easi ly computed. W e refer to a l inear program in the form of equat ion 4.3 as a C O H O linear program. A n d P is cal led a C O H O matrix defined as the ma t r i x w i t h one or two nonzero elements i n each row. 4.1.2 Pr imal -Dual Method C O H O l inear programs can obviously be reduced to s tandard fo rm by add ing ex t ra variables. However, th is w i l l destroy the special s t ructure of or ig ina l C O H O l inear program. W e solve the prob lem using its dua l form that retains the specia l structure. T h e dua l [PS82] of a C O H O L P is a s tandard form l inear p rogram. W e wr i te the C O H O dual linear program or standard form C O H O LP as m in ā bTu s.t. PTu = ETc (4-4) u > 0 and PT has one or two nonzero elements i n each co lumn wh ich is cal led a C O H O dual matrix. 34 T h e p r i m a l / d u a l pair l inear programs have many useful propert ies. B y the Duality Theorem [VanOl , pp.58-67], the p r ima l and dua l l inear programs have the same op t ima l basis and op t ima l value. A n y non-opt imal feasible basis for the dua l prob lem is infeasible for the p r ima l , and produces an upper bound for the p r ima l p rob lem, and vice versa. T h e op t ima l basis is the one that is feasible for bo th p r ima l and dua l problems. Therefore, the C O H O L P can be solved by app ly ing S imp lex to the s tandard form C O H O L P as described in the fol lowing. 4.1.3 Simplex Algorithm There is a c lassical a lgor i thm for l inear programming: S implex . However, it on ly operates on s tandand form l inear programs. Therefore, the C O H O l inear programs are converted into C O H O dua l l inear programs first, then solved by our modi f ied S imp lex solver. In add i t ion , an in i t ia l feasible basis must be prov ied to S implex . F i n d i n g a feasible basis is not t r i v ia l , and it w i l l be dealt w i t h in Sect ion 4.2. S imp lex is a greedy a lgor i thm. D u r i n g each step, it tr ies to replace one co lumn of the current basis w i th a new co lumn to obta in another feasible basis that lowers the cost. L e t B be the current basis, Aq denotes A restr ic ted to th is basis. Le t tj be the co lumn vector defined by: ā¢ tj = A~B1A,J (4.5) where A . j is the j t h co lumn of ma t r i x A . T h e var iable c] = Cj ā Cgtj (4.6) is the relative cost of the j t h co lumn w i th respect to current basis B. If the relative cost is negative, the cost is reduced by in t roduc ing the j t h co lumn into the basis and ev ic t ing some basic co lumn. T h e a lgor i thm must determine which co lumn in the basis to eject when a new co lumn enters. T h e decision is guided by the requirement that the new basis 35 must be feasible and the cost is reduced. T h e index of co lumn to be evicted is T h e act ion of replac ing a co lumn in the o ld basis w i t h a new co lumn to get a more favorable basis is cal led pivoting. S implex pivots to reduce cost un t i l the op t ima l basis is found. T h e cost reduct ion for the pivot is To reduce the cost, costreduce must be negative, therefore ^ should be posit ive. Because P g is feasible basis, wh ich impl ies io,j is non-negat ive, then tf-j. should be negative for a success p ivot ing. If to j is zero, there is no cost reduct ion , and the bases before and after such a pivot are cal led degenerate. Fur thermore , if tjti < 0 , V i = 1, ā¢ ā¢ ā¢ , n , then the l inear program is unbounded, feasible points w i t h arb i t ra ry low cost exist. C O H O l inear program can have degeneracy, but it can not be unbounded because i ts feasible region is an bounded po lyhedron. In the absence of degeneracy, S implex works because at each pivot it reduces the cost whi le there are f inite number of feasible bases. In the presence of degeneracy, S imp lex can loop indef in i te ly through a set of degenerate bases, wh ich is cal led cy-cling. C y c l i n g avoidance is achieved by B land ' s ant icyc l ing a lgor i thm [Bla77],[PS82]. S imp lex can not guarantees f ind ing the shortest pa th f rom the in i t i a l feasible basis to the op t ima l basis, and in fact i t 's an exponent ia l a lgor i thm in the worst case. However, i t per forms quite wel l in pract ice. 4.2 Combination of Simplex and Interval Computation T h e S imp lex a lgor i thm descr ibed above assumes that a l l operat ions are free of errors. Unfor tunate ly , f loat ing point computa t ion introduces round off error. E r ro rs i n the computa t ion may prevent S imp lex f rom produc ing the correct result . A l t hough the k ā arg m m āā (4.7) costreduce (4.8) 36 accumulated error is reduced largely by Laza 's efficient l inear system solver [LazOl , pp.32-36], i t can not be e l iminated. For example, S imp lex depends on comparisons between f loat ing-point num-bers to f ind the new basic co lumn and evicted co lumn. In the presence of i l l -condi t ion ing, the result may be incorrect. C lear ly , a wrong p ivo t ing can make the a lgor i thm fa i l : ā¢ If a wrong basic co lumn is evicted out of the basis, an infeasible basis is reached. ā¢ If a wrong co lumn is brought into the basis, a more cost ly basis is reached that make the a lgor i thm fal l into a cycle of bases. A t the same t ime, C O H O requires an over-approx imat ion result for ver i f i -cat ion. T h e op t ima l value f rom S implex may be under -approx imated because of nearest representat ion round ing off, even though the op t ima l basis is correct. To solve a l l problems above, we combine the S imp lex a lgor i thm and interval computa t ion . A l l values used in the solver are interval numbers rather t han floating point values, as descr ibed i n Sect ion 3.2. T h e interval based S imp lex a lgor i thm returns a set of possible op t ima l bases and a range that contains the op t ima l value. Therefore, over-approx imat ion is guaranteed, and an error b o u n d is prov ided. In add i t ion , the i l l -condi t ion problems are easily detected us ing in terva l ar i thmet ic as descr ibed i n Sect ion 4.4. However, us ing interval ar i thmet ic , the outcome of compar isons of interval numbers can be ambiguous: if the intervals for x and y overlap, then the order ing of x and y cannot be determined. Therefore, the in t roduc t ion of in terval representat ion introduces several problems: 1. P i vo t s might be favourable or not. 2. Bases encountered might be feasible and /o r op t ima l . 37 3. T h e a lgor i thm for f ind ing the in i t ia l feasible basis proposed by L a z a is no longer guaranteed to succeed. T h e fo l lowing describes our interval based solver and a lgor i thm to handle these problems. These methods are not restr icted to C O H O l inear programs, it modif ies the S imp lex a lgor i thm thus works for general l inear programs. 4.2.1 Ambiguous Pivoting and Cycling T h e S imp lex a lgor i thm works by pivot ing to bases that progressively lower the cost unt i l the op t ima l basis is reached. These pivots are determined f rom equations 4.6 and 4.7. However, in terval ar i thmet ic operat ions i n our solver make it possible that there are several ambiguous pivots and no clearly favorable choice. For example, the relative cost may not be clear ly posit ive or negative. T h u s , it might be impossib le to f ind a co lumn to b r ing in such that the new basis is c lear ly favorable. S imi lar ly , it might be impossib le to f ind the basic co lumn to evict. Therefore, c lassical S implex a lgor i thm does not work i f there is no clearly favorable bases because of interval computa t ion . To deal w i t h pivot selection, we uses Laza 's me thod [LazOl , pp.87-90]. T h e a lgor i thm first tr ies to f ind a clearly favourable pivot if possible. If uncerta inty happens dur ing interval comparisons, it tries a l l possibly favorable pivots. In our implementa t ion , there are two steps of a pivot . F i rs t l y , the a lgor i thm finds the co lumn to b r ing i n , determined by relat ive cost f rom equat ion 4.6. It first checks to see i f there is a co lumn whose relat ive cost is c lear ly negative. T h i s involves checking each co lumn that is outside of the current basis un t i l one w i t h clear ly negative relat ive cost is found. E a c h such check requires solv ing one l inear system in l inear t ime, wh ich w i l l be described in the next sect ion. T h u s , the tota l t ime is 0(n(m ā n)) where n is the number of constraints and m is the number of variables. If such co lumn is found, i t is taken, and the other co lumns do not need to be evaluated at this step. If not, then a l l the co lumns w i t h poss ib ly negative 38 relat ive cost are considered. In the second step, the a lgor i thm f ind the basic co lumn to evict , us ing equat ion 4.7. If there is a clear result f rom the argmin funct ion, then it is taken; otherwise, a l l basic columns that might be evicted are considered. A t th is point , we have a set of possib ly favorable pivots, our a lgor i thm branches: a l l are explored, as shown in F igure 4.1. F igu re 4.1: P i v o t Select ion of C O H O L inear P r o g r a m Solver. B ranch ing leads to three types of bases: l ) imfeas ib le basis; 2) op t im ia l basis as expected; 3)more cost ly basis wh ich causes a loop. T h i s branch ing p ivot ing may lead to a more cost ly basis thus make the a l -gor i thm fal l in to a cycle, wh ich precludes the use of s tandard ant i -cyc l ing algo-r i thms [Bla77]. Instead, we main ta in a hashtable of a l l bases that we have v is i ted. If a pivot leads to a basis that we have already seen, it is ignored. It is also possi-ble that infeasible bases are encountered. If a basis is c lear ly infeasible, then it is removed by the a lgor i thm wi thout changing the op t ima l result. Otherwise, it must be qui te close to the feasible region; therefore the error in t roduced is very smal l . B y branch ing, C O H O is guaranteed to vis i t the t ru ly op t ima l basis. T h e branch ing p ivot ing may lead to an exponent ia l a lgor i thm in theory. However, b ranch ing is rare i n pract ice. Thus , the a lgor i thm remains efficient. ā¢ ā¢ ā¢ā=Ā»-J infeasible {^) optimal (~) cycling 39 4.2.2 Possible Optimal Result Simi la r issues arise when test ing bases for opt imal i ty. Ma themat i ca l l y , a basis is op t ima l if and on ly if it is feasible for bo th the s tandard- and C O H O - f o r m L P s . If C O H O encounters a c lear ly op t ima l basis, it is done. Otherwise, each t ime C O H O encounters a basis that is possib ly op t ima l , it checks to see if the corresponding basis for the C O H O - f o r m L P is possib ly feasible. If so, the basis is flagged as possib ly op t ima l . If no c lear ly op t ima l basis is found, C O H O uses the set of possib ly op t ima l bases to determine an interval that it guaranteed to conta in the actua l cost at opt imal i ty . Therefore, we redefine the interface for the l inear p rogram solver. T h e result it returns is a bound for its op t ima l cost, wh ich is usual ly very sma l l , and a series of possible op t ima l bases, wh ich defini tely contains the real op t ima l basis. 4.2.3 Initial Feasible Basis Simp lex requires an in i t i a l feasible basis, whi le f ind ing a feasible basis is non- t r iv ia l for general L P . L a z a [LazOl , 82-89] presented a s imple me thod to f ind a feasible basis. Laza ' s a lgor i thm finds an invert ible basis f irst, then constructs a helper L P : the invert ib le basis is feasible for the helper L P , and the op t ima l basis of the helper L P is feasible for C O H O L P . Laza ' s approach does not work when interval ar i thmet ic is used. A s descr ibed above, the solver can on ly re turn a numer ica l ly possible op t ima l basis for the helper L P thus a possible feasible basis for the or ig inal . C O H O LP. If the basis is actual ly infeasible, then this a lgor i thm fails. Ra the r , we app ly the Big M method [Win87]. For a C O H O Dual LP: min(ciyi + ... + cnyn) yi bm 40 Hi > 0 i = 1 , . . . ,n we add m extra variables zi,..., zm, and make them expensive by assigning this a very large cost, M ; thus the ex t ra variables must be dr iven out of the basis by the no rma l p ivo t ing procedure of S implex. We append an ident i ty ma t r i x to the end of C O H O Dual Matrix A 1: min(ciyi + . . . + Cnyn + Mz\ + . . . + Mzm) s.t. yi an ā¢ din 1 ā¢ 0 yn 0 ā¢ 1 Z\ Q>mn bm Vi>0 i = l , . . . ,n J = 1 , . . , m We can see th is is also, a C O H O Dual LP, it has the same op t ima l value and op t ima l basis as the or ig inal one because the added variables z\,..., zm are expensive and w i l l not appear i n the op t ima l basis. T h i s is equivalent w i t h add ing some redundant constraints, s imi lar w i t h Xi < M, to the C O H O LP. C lear ly , the added extra columns are a clearly feasible basis for th is helper L P . T h e extra column can not appear i n the op t ima l basis. However, there are L P s for wh ich we can not dr ive these columns out no mat ter how big M is. For example, when some 6j = 0, then Zi ā 0, the cost for it is zero and thus can not be dr iven out. So we shou ld not in t roduce such columns dur ing p ivo t ing step; and even if the relative cost of co lumn j is zero, we should t r y to b r ing it into basis, because a l though it can not reduce the cost, it may dr ive out the extra columns w i t h the same cost. A t last, we shou ld remove any poss ib ly op t ima l bases wh ich contain such undesirable co lumns. ^ e r e we assume bi are all nonnegative. If 6, is negative, we just set the ith diagonal element asā1 to make the basis feasible ' ' 41 4.3 Efficient Linear System Solver So far, we have developed an interval based l inear p rogram solver. It is not as effi-cient as the classical S imp lex a lgor i thm, a l though it has several advantages. In the a lgor i thm, so lv ing the l inear systems to compute the tab leau is cr i t ica l for perfor-mance. L a z a [LazOl , pp.32-36] has developed a l inear system solver, tak ing advan-tage of the special s t ructure of C O H O mat r i x , that improves efficiency and reduces accumulated error. We re implemented his a lgor i thm using interval ar i thmet ic and extended the a lgor i thm for bo th C O H O matr ices and C O H O dua l matr ices. We also implemented a hybr id a lgor i thm that improves the accuracy of the result w i t h a smal l t ime penalty. T rad i t i ona l implementat ions of S implex exploi t the proper ty that each basis is a rank-1 update [GV96] of i ts predecessor. T h i s allows the upda ted tableau to be computed in 0{n(n ā m)) t ime where n is the number of constraints in the s tandard form L P , and m is the number of variables. However, th is approach accumulates error f rom one step to the next. Us ing interval ar i thmet ic , the error bounds wou ld qu ick ly become larger than the values themselves. Instead, we explo i t the structure of the P ma t r i x to ob ta in a l inear t ime a lgor i thm for so lv ing the l inear systems that arise in the S imp lex a lgor i thm wi thout error accumulat ion. A s noted earl ier, each row of P corresponds to an edge of a pro ject ion po lygon and thus has either one or two non-zero elements. Accord ing ly , each co lumn of ma t r i x P j has either one or two non-zero elements. N o w , consider the rows of P g . B y the p ivot select ion rules, Pg w i l l never have a row that is a l l zeros. If a row has exact ly one non-zero element, then the value of the corresponding var iable can be determined d i rect ly and the row and co lumn are e l iminated. Such e l iminat ion can be cont inued unt i l every row has at least two non-zero elements. These el iminat ions preserve the proper ty that every co lumn has at most two non-zero elements, and the mat r i x remains square. T h u s , at the end of this phase, every row and every co lumn of the ma t r i x has exact ly two non-zero elements. T h i s ma t r i x has many wonderfu l 42 propert ies; so, we w i l l ca l l th is ma t r i x W. Assume that W is n x n. Consider an undirected graph w i t h vertexes v\...vn such that there is an edge between vertices i and j iff W has a co lumn that has non-zero elements i n rows i and j . B y the structure of W, every vertex i n this graph has degree 2, and the graph can be decomposed into dis jo int , s imple cycles. T h e l inear systems corresponding to these cycles can be solved independent ly. For s impl ic i ty of presentat ion, we w i l l assume that W has a single cycle and note that the general izat ion to mul t ip le cycles is stra ight forward. W e can now permute the rows and columns of W such that Wij 7^ 0 iff j = i or j = {i m o d n) + 1. W e can scale the rows so that every d iagonal element of the ma t r i x is 1. T h e resul t ing mat r i x has the form shown below: W = 1 0 0 -ai 1 0 -a2 0 0 - Ā« n - l 1 (4.9) Now, to solve a l inear system of the form Wu = d, we observe that Ā£ " = i ( < ^ - i ) Ux = Ui+1 = I - E J = I ( Ā„ H ) Pi (4.10) (4.11) where 0k = f l L i ai-F igu re 4.2 summar izes the flow of our l inear system solver. It is straightfor-ward to show that a l l of the steps described above for so lv ing a l inear system can be per formed in 0(n) t ime where n is the number of equal i ty constra ints i n the s tandard fo rm L P . Fur thermore , the on ly step that int roduces cancel lat ion is the denominator 1 ā (3n i n equat ion 4.10. Thus , we can rap id l y ident i fy i l l -condi t ioned l inear systems. 43 Eliminate singleton rows Partition into cycles Figure 4.2: C O H O L inear Sys tem Solver Solve each cycle In the event that a basis is poor ly condi t ioned, we can use equat ion 4.10 to solve for each var iable separately. T h i s results i n an 0(n?) a lgor i thm but avoids accumula t ion of error. In our implementat ion, we use the l inear t ime method first, and if the error exceeds a preset bound , solve again w i t h the more precise 0(n?) algor i thm. In pract ice, the 0(n) a lgor i thm is used almost a l l of the t ime. T h u s we obta in h igh accuracy and low run- t ime. 4.4 Ill-Conditioned Basis and Exception Handling T h e solver descr ibed above works efficiently and re l iably for wel l -condi t ioned C O H O l inear programs. However, when the prob lem is too bad ly condi t ioned, for example, when Pn is close to 1 in equat ion 4.10, the error bound of the l inear system solver might be so large that the op t ima l value is overapprox imated too much to be useful. In these cases, a s l ight ly pr imal- infeasible but wel l -condi t ioned basis is a better overapprox imat ion of the op t ima l cost. 4.4.1 Conditioning Interval computa t ion provides and error bound that makes it qui te easy to detect bad ly condi t ioned problems. However, this method requires solv ing the l inear sys-tems even though it is bad ly condi t ioned. We need a fast me thod to detect bad ly 44 i i l l -condi t ioned problems. T h e cond i t ion number of a mat r i x is the rat io of the largest to smallest singular value i n the singular value decomposi t ion, it is an est imate of how many digi ts are lost in so lv ing a l inear system w i t h that mat r i x . B u t it is also expensive to compute the exact cond i t ion number. A n approx imat ion technique is required. L a z a proposed an a lgor i thm to compute the upper b o u n d on the condi t ion number of C O H O ma t r i x [LazOl , pp.41-46]. However, the est imat ion is often much larger t han the actua l condi t ion number. Thus , problems can be needlessly identi f ied as bad ly condi t ioned at beginning. Fur thermore, the est imat ion for a ma t r i x and its t ranspose can be different, wh ich makes it possible that a basis is wel l -condi t ioned for C O H O dua l l inear program whi le i l l -condi t ioned for C O H O l inear program. In the imp lementa t ion of solver, unexpected except ion might be th rown at the f inal step that checks the feasibi l i ty of bases for C O H O l inear programs. F ina l l y , the computa t ion of cond i t ion number should be much faster t han the l inear system solver, otherwise approx imat ion is neither efficient nor necessary because interval ar i thmet ic provides error bound . Laza ' s computa t ion is not easy enough compared w i th his fast l inear system solver. We developed a new a lgor i thm for comput ing lower bound of the condi t ion number of a C O H O mat r ix . It catches most bad ly condi t ioned problems. W h e n a bad ly condi t ioned prob lem skips through, the intervals of the result are very large, and the prob lem is then detected. We use the rat io of the largest to smallest eigenvalue of the mat r i x as the underapprox imat ion . T h e eigenvalue for C O H O ma t r i x or C O H O dua l ma t r i x is easy to compute. Le t A denotes the eigenvalues of C O H O ma t r i x or C O H O dua l ma t r i x A , we have ( 1 - A ) n = n a i (4.12) i=i therefore, the eigenvalues are A = l - / ? y n ( c o s ā + i - s i n ā ) (4.13) n n n 45 and our est imat ion is max( |A i | ) = i+p1/n T ^ " zā it n is even \i-Pn/n\ ^ / l + V " - 2 P ^ n c o s 2 Ā£ l i Z 2 1 . ( 4 - 1 4 ) _i ii I T n I Q P V P T I min(|Aj | ) N o w , let us prove it is the lower bound of cond i t ion number . TT- i f is even \l-Pn' I Theorem 4 .4 .1 F o r a nonsingular matrix A, let A i , ā¢ ā¢ ā¢ , A ā be the eigenvalues of A and o~i, ā¢ ā¢ ā¢, an be the singular values, then < minfg,'^ Proof Suppose A is one eigenvalue of A, then 3 an eigenvector x ^ 0 such that A x = Ax . B y the proper ty of no rm A| |x | | = | |A| | . | |x | | = | | ^ | | < | | y l | | . | | x | | (4.15) therefore, A < P|| (4.16) and it is st ra ight forward that max( |Ai | ) < PH (4.17) S imi lar ly , for ma t r i x A"1, have m a x d A - 1 ! ) < p-1!! (4.18) and because eigenvalues of A^1 are reciprocal of that of A, thus PH.p-1!! > m a x ( | A i | ) . m a x ( | A r 1 | ) mm(|Aj | ) B y def in i t ion max^ ) = pH (4.20) and ; min ( (Ti) = m a x f f f - 1 ) = p- J|| (4.21) 46 Therefore, max(|Aj | ) < m a x ( a i ) min(|Aj|) ~ min(<7;) I 4.4.2 Use of Non-optimal Bases Once a bad ly condi t ioned l inear system detected, we have to f ind an good over-approx imat ion for the i l l -condi t ioned op t ima l basis. A n good approx imat ion must satisfy: 1) the approx imated cost is close to the op t ima l cost, 2) the corresponding basis is wel l cond i t ioned 3) it is over approx imated. < A square ma t r i x is i l l condi t ioned if at least one row (column) is nearly a l inear combinat ion of the other rows (columns). T h e number of rows (columns) that are near ly equal to l inear combinat ion of the remain ing rows (columns) represents the degree of i l l -condi t ion. Note that changing any co lumn in the i l l -condi t ioned cycle changes the value of (3n and can remedy the prob lem. T h i s s i tuat ion is quite rare; so, we s imp ly d iscard one of the columns f rom the s tandard form l inear program and the l inear p rogram solver is restarted. T h i s may result i n a sub-op t ima l so lut ion to the s tandard form L P and therefore a super-opt imal so lut ion to C O H O - L P . T h i s preserves our guarantee of overapproximat ion. To avoid excessive growth i n the reachable region, we test the discarded constraint at the end of the L P solut ion. If i t is c lear ly v io la ted, then we br ing that co lumn back into the L P and reject a different one. W e cont inue i n this fashion unt i l we f ind a good qua l i ty op t imum. T h e method has a geometric interpretat ion. For C O H O l inear programs, basic co lumns (rows) i n the dua l (pr imal) represent inward halfspace normals in the p r ima l . T h e feasibi l i ty of a basis i n the dua l means that the p r ima l cost vector is a posit ive combinat ion of basis columns. Therefore, the p r i m a l cost vector lies inside the cone generated by the basic inward halfspace normals , also cal led the basic cost cone. If the i l l -condi t ioned basis is non-opt imal , then our me thod w i l l not affect 47 the f inal result . Le t (3 be the set of feasible bases of the C O H O dua l l inear program. T h e n once a bad ly condi t ioned basis found, some basic co lumn c is exclude f rom the feasible basis. Therefore, the set of feasible bases 0' for,the modi f ied l inear program is the subset of (3. T h u s , our method on ly reduce the feasible region and the number of feasible bases. B y the assumpt ion that the i l l -condi t ioned basis is not op t ima l , there exists a basic co lumn c that is not i n the op t ima l basis. Therefore, the op t ima l basis is not changed after exc lud ing co lumn c. O n l y an i l l -condi t ioned op t ima l basis matters. N o w consider a two-dimensional C O H O l inear program. Here, there are only two basic co lumns, i l l -condi t ion ing indicates the two normals are near ly paral le l , wh ich means the op t ima l vertex is either h ighly obtuse or h igh ly acute, as shown in F igure 4.3. Because a l l bases v is i ted in the S implex a lgor i thm are feasible for the d u a l 2 , the op t im iza t ion d i rect ion of the p r ima l lies inside the op t ima l cost cone. Therefore, when the op t ima l vertex is very obtuse, i.e. the op t ima l cost cone is very acute, then the op t im iza t ion d i rect ion is near ly perpendicu lar to sides of the op t ima l vertex. T h i s suggests that vert ices that l ie on the l ines cor responding to the two constraints are good choices for approx imat ing the opt ima l cost. One constraint in the basis is d ropped because it is redundant . T h e approx imated vertex is obta ined by replace this constraint w i t h some other constraint out of the basis in our a lgor i thm. For the h igh ly acute vertex case, it is precluded in the geometry phase of C O H O , wh ich replaces such polygons w i th a smal l rectangle. T h u s i t is not consid-ered. N o w , let us t r y to extend the idea to higher d imensions. C o n s i d e r a C O H O l inear p rogram of d imens ion d whose op t ima l basis B is i l l -condi t ioned. Let us assume the degree of i l l -condi t ion is one, it is 's t ra ight forward to generalize to higher 2 For our solver, it is possible that little infeasible is visited. However, if the infeasible basis ill-conditioned, it is safe to drop it because the basis is infeasible which should not be visited at all. 48 Figure 4.3: Types of O p t i m a l 2D Vert ices: l ) h i gh l y obtuse op t ima l vertex and h ighly acute cost cone; 2) h ighly acute opt ima l vertex and h igh ly obtuse cost cone. degrees. Le t Gk be the l inear subspace generated by a l l basic co lumns except the fcth: {B\, ā¢ ā¢ ā¢, Bk-i, Bk+\, ā¢ ā¢ ā¢ Bd}- T h e subspace is d ā 1 d imens iona l , i.e., a hyperplane. E a c h of these hyperplanes generated in the manner is a face of the op t ima l cost cone. Assume that the last basic vector Bd is near ly a l inear combina t ion of the other co lumns, whi le the others are independent. Therefore, the pro ject ion of Bd onto direct ion rid, wh ich is the no rm of Gd and also the intersect ion of faces corresponding to the 1 s t , ā¢ ā¢ ā¢ , d ā Ith constraints, is almost zero. A s shown in F igu re 4.4, if the basis vector Bd onto the hyperp lane Gd is posi t ive, then the cost cone is h igh ly acute and near ly d ā 1 d imensional . Therefore, the op t im iza t ion d i rect ion v is almost on the d ā 1 d imensional hyperp lane; it is almost perpendicu lar to vector rid- Therefore, any vertex on the no rm of Gd is a good approx imat ion . T h e dth constraint is dropped to f ind the approx imated vertex s imi lar w i t h the two d imensional case. Otherwise, i f the project ion of basis vector Bd onto the hyperp lane Gd is negative, the cost cone is h igh ly obtuse and it is also near ly d ā 1 d imens iona l . T h e angle between Bd and the no rm of Gd is h ighly acute. However, i n pract ice such 49 F igu re 4.4: H igh ly Obtuse O p t i m a l 3D Vert ices cases are removed s imi la r ly in the geometry step. Therefore, our method is reasonable and prac t ica l . T h e error in t roduced by the approx imat ion is qui te smal l i n our appl icat ions. 4.5 Implementation T h e C O H O solver is implemented in coho.lp and coho.Ip.solver packages. T h e Ip pack-age defines the interface for a l inear program, inc lud ing classes for l inear program, constraints, basis and result . It also defines a interface for l inear p rogram solver, each solver should implement it. T h i s architecture makes it qui te easy to change to a new l inear program solver. T h e solver package implement the l inear system solver i n the CohoMatr ix class and the l inear program solver i n the CohoSolver class. T h e solver is integrated w i t h interval number implemented i n coho.common.number package descr ibe i n Chap te r 3.2. T h e CohoMatr ix is a subclass of Doublelnterval-Matr ix i n coho.common.matr ix package. T h e st ructure of C O H O solver is shown i n F igu re 4.5. 50 Ā« L P Ā» LP Constraint LPbasis LPresult LPsolver LPError Ā« I n t e r f a c e Ā» LPSolver D-1p: LP +0Pt() LP -eq: Constraint -neq: Constraint -costVector : Matrix +dualLP<) +createCoho(j +createCohoDua\() Ā« S o " L v e r Ā» CohoMatrix CohoSolver Ā« M a t r i x Ā» DoublelntervalMatrixl' CohoSolver ā¢Ip: LP +opt() - CohoMatrix I +linearSystemSolver() 1 I Ā«Numbe r Ā» CohoDoublelntervall Figure 4.5: C lass Hiearchy of C O H O L inear P r o g r a m Solver In sect ion 4.2.3, the Big M method is descr ibed. T h i s a lgor i thm guarantees f ind ing an in i t ia l feasible basis for any l inear programs, however, i t does not use the special s t ructure of C O H O l inear programs. In the implementa t ion , a fast method is appl ied f irst, the B i g M method is used on ly if it fai ls. T h e idea is very s imple: For the s tandard C O H O l inear program, there are at most two nonzero elements in each co lumn, and usual ly there are many co lumns that have on ly one nonzero element. Therefore, the solver tr ies to f ind a basis that there is on ly one nonzero element i n each row and each co lumn. It 's obviously that th is k i n d of bases is invert ib le and feasible. If such a feasible basis is found, no ex t ra variables need to be in t roduced, wh ich makes the fol lowing p ivot ing much easier and faster. If Big M method is per formed, ex t ra variables are added to the l inear pro-51 gram. There might be several i l l -condi t ioned bases encountered du r ing the p ivot ing to dr ive out these ex t ra variables. Because the weight for ex t ra variables are expen-sive, the in i t ia l feasible basis may lead to some bases that have very large cost whi le s t i l l c lear ly favorable for p ivot ing, for example, two near ly para l le l edges for the 2D case. However, th is w i l l in t roduce some unexpected i l l -condi t ion bases, wh ich slows our a lgor i thm. To avoid such k ind of bad ly condi t ioned basis, ex t ra variables are never al lowed to enter into the basis once dr iven out in the p ivot ing. A n d if a new favorable basis is obta ined by d r i v i ng out one ex t ra var iable, the cond i t ion ing of the new basis is checked. If i t is h ighly i l l -condi t ioned, the a lgor i thm tries to f ind a new wel l -condi t ioned whi le clearly favorable basis as possible. T h i s removes lots of redundant bad ly condi t ioned bases. T h e C O H O solver is complex to debug, it depends on several packages, inc lud-ing the number and matrix packages. However, it is easy to check the opt imiza t ion of a basis. T h e solver and related packages are val idated by checking the feasibi l i ty of bases re turned by solver for bo th p r ima l and dua l l inear programs. 52 Chapter 5 Projection A t each step of C o h o , the advanced projectahedron must be projected onto the 2D subspace, as descr ibed in Chap te r 2. E a c h t ime advanced slab is descr ibed as a series of l inear inequal i t ies, i.e. the feasible region of a C O H O l inear p rogram. P r o d u c i n g the pro ject ion polygon f rom these l inear inequal i t ies requires solv ing a series of l inear programs using our C O H O solver. Coho ' s a lgor i thm for comput ing project ion polygons is in t roduced in sec-t ion 5.1. Sect ion 5.2 describes the issues we encountered for a prac t ica l implementa-t ion and reduc ing error. O u r a lgor i thm for convert ing a general pro ject ion prob lem to the s imple one is descr ibed in section 5.3. Sect ion 5.4 analyzes the error of the C O H O pro ject ion a lgor i thm. 5.1 Algorithm A t the end of a t ime step, a slab is advanced to a pro jectahedron according to the l inear mode l . To main ta in the projectahedron representat ion, the project ion onto the or ig ina l 2 D subspace is needed. Formal ly , the prob lem is to project a pro jectahedron descr ibed as: P-x>b (5.1) 53 where P is a C O H O ma t r i x , onto project ion plane V, generated by basis vectors e j ,e j Ā£ {||ej|| = ||e,-|| = l , e , ā¢ ej = 0}. It is a special case of comput ing the project ion po lygon of a po lyhedron onto a 2D subspace. T h e basic idea is to f ind points on the pro ject ion po lygon by max im iz ing the value along each di rect ion on the project ion plane. T h e op t im iza t ion prob lem is to solve a l inear p rogram, of wh ich the feasible region is the pro jectahedron and the op t im iza t ion d i rect ion is the vector on the pro ject ion p lane, for example, the norma l of an edge. It can be solved by the C O H O solver. However, there are an inf in i te number of possible direct ions on the pro ject ion plane. It is impossib le and unnecessary to solve the opt imiza t ion problems for every possible d i rect ion. To reconstruct a po lygon, on ly the polygon vertices are required. A s shown i n F igu re 5.1, so lv ing the opt im iza t ion problems along the normals of po lygon edges is enough to f ind a l l po lygon vertices and the project ion po lygon. t F igure 5.1: P ro jec t i on A l g o r i t h m Us ing L inear P rog ramming : po lygon vertices are found by op t im iz ing along the normal of edge. O u r a lgor i thm finds the normal for each edge of the pro ject ion polygon in a counter clockwise sweep. G i v e n the right endpoint of the previous edge, i.e, the left endpoint of current edge, the normal of the current edge is computed using the feasibi l i ty of op t ima l basis returned f rom previous l inear p rogram; then a l inear p rogram, of wh ich the opt im iza t ion d i rect ion i s t h e no rma l , is const ructed, the right 54 endpoint of current edge is found by solving it. T h i s step, as shown i n F igu re 5.2, is per formed unt i l a l l edges are found, then the pro ject ion po lygon is reconstructed easily. F igu re 5.2: F i n d the N o r m a l of Cur rent Edge Us ing O p t i m a l Bas is of Prev ious L inear P r o g r a m T h e essential part of our a lgor i thm is the currNorm func t ion wh ich computes the no rma l of the current edge. G i ven the l inear p rogram, of wh i ch opt imiza t ion d i rect ion is the no rma l of the previous edge, the op t ima l basis B that corresponds to the left endpoint of the current edge is op t ima l for any l inear programs w i th op t im iza t ion d i rect ion between the previous normal and the current normal. A s shown i n F igu re 5.2, when the opt imiza t ion d i rect ion exceeds the current normal, B is no longer op t ima l . It becomes a feasible but non-op t ima l basis for the C O H O l inear p rogram and thus infeasible for the C O H O dua l l inear p rogram. Therefore, the current normal is the cr i t ica l d irect ion that forces B to be infeasible for the s tandard form C O H O l inear program. Fo r the s tandard form C O H O L P , the infeasible basis forces at least one var iable to become negative. Le t e, + aej be the op t im iza t ion d i rect ion, for a basis B, the variables x are given by: x = PBT-(ei + aej) = P^ei + aP^ej (5.2) = 7r + an 55 where TT = PB r e , and n = aPB Tej. Obviously , the a value of current normal is the smallest that makes some var iable Xi of x nonposi t ive: CtcurrNorm ~ Ulin( ) (5.3) Ā» Vi Here, 7T^ shou ld be posi t ive and rji should be negative. Once the d i rect ion exceeds current normal, the a value is greater than acurrNorm, then the var iable Xi is negative , and B is no longer feasible. O f course, any op t im iza t ion direct ion which is ahead of the previous normal also makes the basis B infeasible. T h i s can be easi ly avoided by forcing the a value to be greater than the otcurrNorm. G i v e n an op t im iza t ion d i rect ion, the corresponding l inear program is solved; bo th the op t ima l basis and op t ima l value are returned by the C O H O solver. A po lygon can be represented as a col lect ion of edges or a set of vert ices. Here, the first me thod is used and an edge is represented by the op t ima l d i rect ion and its op t ima l value. T h e vertices are reconstructed f rom edges at the final step; rather comput ing f rom op t ima l bases direct ly, a l though they are usual ly the same. There are several reasons for choosing this approach. F i r s t , the op t ima l basis returned f rom the C O H O solver is possibly op t ima l , it might not be the exact vertex, even though it w i l l be qui te close. Second, if some edge is omi t ted , our me thod guarantees the final result is over approx imated , whi le it is under approx imated us ing the vertices direct ly. A s shown in F igure 5.3, if the middle edge is sk ipped , our method w i l l replace the 1 s t and 2 n d vertices w i t h the 3 r d vertex. Otherwise, the 2 n d vertex is ignored and an unacceptable under approx imat ion is p roduced. T h e basis corresponding to r ight endpoint of current edge should be found to compute the no rma l of the next edge, as shown in F igure 5.2. Ideally, if the op t im iza t ion d i rect ion is over the norma l , the op t ima l basis for the l inear program is the one corresponding to the r ight endpoint . However, for our C O H O solver, only possibly op t ima l bases are available. Therefore, to dr ive out the left endpoint f rom the set of possibly op t ima l bases, the opt imiza t ion d i rect ion shou ld be sl ight ly 56 under approximation when the vertex of the middle edge is skipped 1 2 3 over approximation when the normal of the middle edge is skipped F igu re 5.3: Overes t imat ion of Pro jec t ion Po lygon W h e n an Edge is Sk ipped greater than the current normal. O f course, th is may sk ip some edges, but it only introduces sma l l amount of over approx imat ion, wh ich w i l l be proved in section 5.4. 5.2 Implementation T h e pro ject ion plane is par t i t ioned into four quadrants by basis vectors ej and ej. T h e norma l of po lygon edge is represented as combinat ion of basis vectors. A s shown in F igu re 5.4, i t 's ej + aej in the first quadrant , ej ā aej i n the second quadrant and so on , and a is always posi t ive. T h e reason for this representat ion is that i t 's easy to compute the no rma l of the edge, as shown i n equat ion 5.3. F igu re 5.4: T h e Representat ion of Op t im i za t i on D i rec t ion 57 Because of round off error and error f rom the C O H O l inear p rogram solver, the op t im iza t ion d i rect ion might not be the exact no rma l of the po lygon edge. If the adjacent normals are near ly the same, i.e. the adjacent po lygon edges are nearly on the same l ine; i t 's better to replace the two normals w i t h one and combine the two edges as one. T h i s w i l l in t roduce smal l error, wh ich w i l l be examined in section 5.4. In the imp lementa t ion , we define a threshold e, the currNorm funct ion w i l l force the angle of normals of current and previous edges to be bigger than a smal l amount e. If the angle is less than e, these two normals are combined and one is d ropped. W h e n the norma l is close to the second basis vector, i.e. a is close to inf ini ty, it is replaced w i t h the second basis vector, wh ich prevents large computa t ion error caused by a huge value for a. A s in i t ia l i za t ion of our a lgor i thm, the first op t im iza t ion d i rect ion is required. A l inear p rogram w i t h opt imiza t ion direct ion ej is solved first and the basis for the left endpoint of the first edge is found. Us ing this basis, the no rma l for the first edge is found, as shown in F igure 5.5. F igu re 5.5: F i n d the In i t ia l Op t im iza t i on D i rec t ion for P ro jec t i on A l g o r i t h m However, the in i t ia l op t im iza t ion direct ion might be the exact no rma l of the 58 first edge, as shown in F igure 5.5. In this case, the first edge is omi t ted . Therefore, at the end a lgor i thm enters into the first quadrant again and check whether i t 's an exact no rma l or not. If i t is, the first edge is added. To f ind the right endpoint , the normal of current edge is moved forward by a smal l amount e, then a l inear program is solved to max imize value along this d i rect ion. T h i s may omi t some edges, but the error is smal l . However, the left endpoint may also be in the set of possible op t ima l bases. We w i l l t u rn the direct ion increasingly unt i l the left endpoint is dr iven out or the change in angle exceeds a m a x i m u m angle. If the no rma l of the current edge is i n the next quadrant , increasing a w i l l not change the op t ima l basis. T h u s the currNorm can not f ind a value for a that makes the op t ima l basis infeasible for the standard C O H O l inear p rogram, the program should j u m p to the next quadrant and reset the value a. T h e in i t ia l a should be zero; however, in the implementat ion, it is assigned w i t h a negative value. T h e currNorm funct ion forces a to be non-negative. However, the nextNorm funct ion might re turn ' jump to next quadrant ' incor-rectly. T h i s is caused by the C O H O solver. A s descr ibed above, the r ight endpoint of the current edge is required to compute the norma l of the next edge. However, the left endpoint might be i n the set of possible op t ima l bases. If the left endpoint is used, the a value for current norma l is returned again, i t is not bigger than the current a by the requi red amount , e, therefore, the a lgor i thm j u m p s to the next quadrant . T h i s specia l case w i l l in t roduce large amount of over approx imat ion error. To reduce i t , back t rack ing is employed when j ump ing to the next quadrant . A s shown i n F igure 5.6, the angle between current norma l and the second basis vector is bisected, the new di rect ion is used as the opt imizat ion d i rect ion and the op t ima l vertex is found. If the currNorm per formed correct ly and there is no edge omi t ted , the result should be the same w i t h that of the previous l inear p rogram, then a lgor i thm jumps 59 jump to next quadrant i : Figure 5.6: Back t rack ing to Reduce E r ro r W h e n J u m p to Nex t Quadran t to the next quadrant . Otherwise, some edges were sk ipped and a new op t ima l basis that is far f rom the current op t ima l basis should be found. T h e a lgor i thm tries to f ind the edges omi t ted : it bisects the angle between current no rma l and the second basis vector and moves the opt imiza t ion d i rect ion back toward the current norma l un t i l the l inear program has the same op t ima l basis w i t h a current basis or the angle is not bigger t han angle of current norma l by a smal l amount e. There is also a special case when the a lgor i thm j umps two quadrants in a single step, for example, the normal of current edge is in the first quadrant whi le the no rma l of next edge is in the th i rd quadrant . In such case, currNorm might per form incorrect ly . T h e solut ion is to preclude th is specia l case. A s shown in F igu re 5.7, an ex t ra edge in the middle quadrant ax + by = c is added. T h e normal of the edge is y ā x wh ich is obviously in the middle quadrant and the value of c is v0pt ā vConst where vconst is a proper posi t ive constant and vopt is the op t ima l value for the l inear programs w i t h op t ima l d i rect ion as y ā x. A s the result, there are 60 no two quadrants j umps and the feasible region is reduced. However, th is is not an under approx imat ion because at the reconstruct ion step as describe i n the fol lowing paragraph, th is edge is removed. F igu re 5.7: Spec ia l Case of J u m p i n g T w o Quadran ts Con t inuous ly A t last, a col lect ion of edges is found. To get the vert ices of the po lygon, the intersect ion point of adjacent edges are computed. T h e y are solved by our l inear system solver efficiently because they are C O H O l inear systems. 5.3 Converting from End-of-Step to Beginning-of-Step Projection T h e previous sect ion descr ibed project ion for pro jec tahedra descr ibed by equa-t ion 5.1. A f te r mov ing forward according to the l inear mode l add an edgi -3Ā»-X x' = Mx + q (5.4) 61 the forward pro jectahedron is describe as P-E~lx>b (5.5) where E = e M A t is the forward Eu le r operator. In general, PE~l is no longer a C O H O ma t r i x wh ich prevent us f rom using our fast l inear system solver to compute 7r, 77 and reconstruct the vert ices at the last step of pro ject ion. T h e solut ion is convert ing the problem of pro ject ing P ā¢ E~lx > b onto the pro ject ion p lane generated by basis vectors ej and ej to the prob lem of pro ject ing Px > b onto pro ject ion plane generated by basis vectors E e , and Eej. T h e poly-gon vertex for the or ig inal p rob lem p and its corresponding po lygon vertex for the converted prob lem p' have the relat ion p = Ep (5.6) T h e me thod has an algebraic exp lanat ion proposed by M a r i u s [LazOl]. In the pro ject ion step, given an opt imizat ion direct ion d, a l inear p rob lem m i n ( d r ā¢ x) s.t. P-E'lx>b should be solved. Le t y = E - 1 ā¢ x we have m i n ( ( Ā£ T ā¢ d)T ā¢ y) s.t. P-y>b Therefore, l inear programs for or ig inal and converted pro ject ion problems have the same op t ima l value and yopt = E~1x. O u r method is va l id . T h e method also has a geometric exp lanat ion. Le t us assume that a polyhe-dron is represented as P ā¢ E~lx > b i n the coordinate C w i t h basis (e i , ā¢ ā¢ ā¢ , en)T, and the pro ject ion plane is generated by ej and ej. Now, i n the new coordinate C w i th new basis defined as E(e\, ā¢ ā¢ ā¢, e ā ) T , the representat ion for the same projectahedron 62 should be Py > b, and the basis vectors for the pro ject ion plane should be Eef and EeJ. In the new coordinate C, the vertex of pro ject ion po lygon is computed as p', its representat ion i n coordinate C should be p = Ep'. Therefore the two problems must have the same project ion po lygon because they are s imp ly different representat ions of the same po lyhedron and project ion plane. F igure 5.8: Convers ion between End-of -Step Pro jec t ion and Beginning-of -Step P r o -ject ion F igu re 5.8 i l lustrates this conversion. In the example, the basis for the coor-d inate system C is and (ei,e2,e 3) 1 0 0 \ 0 1 0 V Ā° 0 l ) (5.7) / 1 0 0 - 1 0 0 0 1 0 V 0 - 1 0 / ,E 2 0 0 \ 0 1 0 y o 0 1 / ,6 = / - l \ - 1 - 1 v - 1 / (5.8) T h e pro ject ion plane is the one generated by ei,e2- T h e n the new coordinate C should be: (ei,e 2,e' 3) = E(ex, e 2 v e 3 ) 2 2 0 0 \ 0 1 0 V 0 0 1 / (5.9) 63 T h e po lyhedron is of the simple form P > b. T h e basis for pro ject ion plane in C is changed to e[ = ( 2 , 0 , 0 ) r and e'2 = ( 0 , 1 , 0 ) T . N o w , i n C it is easy to see the project ion po lygon w i t h vertices / 1 - 1 - 1 1 \ 1 1 - 1 - 1 (5.10) y o o o o J N o w mu l t i p l y i t w i t h E to obta in p= Ep' = { 2 - 2 - 2 2 \ (5.11) 1 1 - 1 - 1 \ 0 0 0 0 J which are exact ly the vertices of project ion polygon i n the or ig inal coordinate C. 5.4 Error Analysis O u r pro ject ion a lgor i thm overestimates the project ion po lygon , however, the error in t roduced is qui te smal l . F i r s t , the C O H O l inear system solver returns possib ly op t ima l vertices. Such a vertex might be infeasible for the C O H O l inear program. Assume the m a x i m u m distance f rom the poss ib ly op t ima l vertex to the feasible region of C O H O l inear program is 7. In other words, the C O H O solver can clear ly detect the infeasibi l i ty of points for wh ich the distance to the feasible region is greater t han 7. A s shown i n F igu re 5.9, the project ion a lgor i thm finds a vertex f irst, then the corresponding edge is computed ; w i th the normal of the edge, a new l inear program is constructed and i ts op t ima l vertex is solved as the vertex i n the next step, and so on. T h e f inal po lygon computed might be different f rom the correct one. It might have some edges that are close to the correct instead, or have some ex t ra edges. However, if no edge is sk ipped, a l l po lygon vertices of the overest imated polygon are op t ima l vert ices of corresponding l inear programs, and the distance to the feasible 64 region is less than 7. It is obvious that the m a x i m u m increase in area caused by the C O H O solver is 7 L where L is the perimeter of the pro ject ion po lygon. v V feasible region possible feasible region projection polygon vertex of polygon normal F igu re 5.9: P ro jec t i on E r r o r Int roduced by C O H O L inea r P r o g r a m Solver However, i t is possible that some edges are sk ipped. In th is case, the polygon vertex might be new point created by intersect ing two edges, rather than op t ima l vertex of C O H O l inear program. Such vertices are c lear ly infeasible; i n wh ich cases the error is not bounded by jL as described above. There are three cases in which edges might be sk ipped: 1. A s desired, our a lgor i thm skips the edge if the angle between it and the previ-ous edge is less than e, as shown i n F igu re 5.10. It is easy to compute that the m a x i m u m possible error in t roduced is ^ 2 s i n e Ā« ^ 2 e , where I is the length of edge sk ipped. Therefore, the m a x i m u m to ta l error of th is k i n d is less than 2. Not ice that our a lgor i thm uses the left endpoint to f ind the no rma l of the edge. However, the C O H O solver might return points close to the left endpoint as the op t ima l vertex, and the normal computed might differ f rom the correct one. If 65 feasible region projection polygon over approximated area d normal previous edge ^ F igu re 5.10: P ro jec t ion E r r o r Int roduced by Ignor ing Sma l l Ang le it is smal ler , no edge is sk ipped, but some redundant edges may be in t roduced. If i t is greater, some edges might be sk ipped. However, the computed normal is greater than the correct no rma l by no more t han (3 = arc t an ^ as shown in F igu re 5.11. U s u a l l y th is angle is smaller than e, and the error is inc luded in U 2 s in e ā¢ descr ibed above. If the angle is greater than e, then the edge is ' ' ' ā¢ ā¢ ā¢ U r r r T i T . - ^ - o v e r a p p r o x i m a t e d area F igu re 5.11: P ro jec t ion E r ro r Introduced by Sk ipp ing Short E d g e very short because f rom ^ Ā« /3 = arc t an ^ > e, we-know V is bounded by 2 wh ich is a quite smal l value in pract ice. It can be computed that the tota l error of this k i n d is bounded by \^L. In fact, the error is counted i n the error fL caused by the C O H O solver; because no new points are in t roduced in the computed po lygon i n this case, a l though some edges are sk ipped. A l l vertices are f rom the C O H O solver 's op t ima l result. Therefore no vertex is far f rom the feasible region by 7, each vertex is in the possible feasible region of C O H O solver. e normal feasible region possible feasible region projection polygon 66 3. A l though t i t occurs rarely, the C O H O solver might re tu rn the left endpoint (or nearby) of the previous edge as the op t ima l vertex. In this case, the norma l computed might be less than the previous no rma l , hence a lgor i thm tries to j u m p to the next quadrant . A s described earl ier, back t rack ing is employed to reduce the error. It stops on ly if there is new vertex found if the angle is greater t han the previous normal by 2 a , whi le there is no new vertex if it is greater by on ly a. In this case, the m a x i m u m possible area increase is less than ^l2 s in 2 a , where I is the distance between the new vertex and previous vertex, as shown i n F igure 5.12. Fur thermore , I is less than a smal l constant c, otherwise, bisect ion is performed between a and 2 a to f ind the edges ommi t ted . Therefore, the area increase by j u m p i n g to next quadrant is bounded by a smal l constant. In fact, if efficiency is not considered, bisect ion can be per formed unt i l no edge is sk ipped, then no error is in t roduced. ā feasible region - - projection polygon -> bisection direction ^ overapproximated area no new vertex found, |V I v I v a Figure 5.12: P ro jec t i on E r r o r Int roduced by J u m p i n g to Nex t Quadran t Therefore, the m a x i m u m possible error for our pro ject ion me thod is about 7-L + \ l? s in e. e is a constant in the algor ihtm and it is set to l e ā12 in the current imp lementa t ion . W e believe 7 is quite smal l based on the current exper imenta l 67 new vertex found result repor ted i n Chap te r 6. T h u s the error is quite smal l i n pract ice. 68 Chapter 6 Experiments T h i s chapter shows two examples of C O H O : a two d imensional S ink example and a three d imens iona l V a n der P o l osci l lator; these demonstrate the correctness of C O H O and robustness of the implementat ion. T h e performance of the l inear p rogram solver and pro ject ion a lgor i thm is analyzed f rom the exper imenta l da ta . 6.1 Application of Coho Verification System W e first appl ied C O H O to the Sink example f rom [DM98] , a two-d imensional , l inear system. T h e example has the same dynamics mode l w i t h D a n g and Ma le r ' s : x = ā2x ā 3y (6-1) y = 3z - 2y However, to prove that C O H O has the capabi l i ty to deal w i t h more complex system, we have changed the in i t ia l region f rom the rectangle w i t h the d iagonal [(0.1,0.1), (0.3,0.3)] . to a po lygon w i t h vertices [(0.1,0.1), (0.3,0.1), (0.4,0.4), (0.1,0.3)]. (6.2) T h e reachable space is shown i n F igure 6.1. T h e result is much more accurate compared w i t h the result f rom F igure 5.4) in [DM98] . O u r result c learly shows the boundary and d is t inct cycles of the 69 Figure 6 . 1 : T h e Resu l t of a T w o Dimens iona l L inear M o d e l E x a m p l e : S ink sp i ra l , whereas D a n g and Ma le r ' s result merges them together and has much more overapprox imated space. T h e result is even better t han that of the earl ier version C O H O , as shown in F igu re 4 of [GM99] . In the earl ier vers ion, the sp i ra l converges to a smal l c i rcu lar shape; whereas in the current vers ion, the reachable set of the last step reaches a much smal ler po lygon that is near ly a point as seen in the F igure 6 . 1 . T h i s demonstrates that the over approx imat ion is much smal ler than before. T h e greater accuracy of new C O H O arises ma in l y f rom two factors. F i r s t , the C O H O l inear program solver offers t ighter over approx imat ions of the op t ima l value, wh ich is ma in l y because of i ts hand l ing of bad ly condi t ioned problems. Second, the pro ject ion a lgor i thm is more accurate; it skips fewer edges t han the earlier version and thereby produces more accurate project ion polygons. A l t h o u g h the S ink example is two-d imensional and thus no project ion is tak ing place; the project ion 70 a lgor i thm is per formed to compute the polygon representat ion of the reachable space f rom the l inear program wh ich represents the advanced face at step 4(a)iv descr ibed i n Sect ion 2.4. T h e reachable space of the S ink example above is two d imens iona l and the system dynamics are l inear. Therefore, it does not use the projectagon technique. We now app ly C O H O to a three-dimensional , non- l inear mode l to demonstrate the reachabi l i ty ca lcu la t ion using projectagon techniques. T h e mode l is desired f rom a V a n der P o l osci l lator ( 3 V d P ) [HS74, pp. 217f] where we have added the z var iable to provide an example w i t h three dimensions and demonstrated C O H O ' S project ion capabi l i ty. T h e mode l is descr ibed by the fol lowing nonl inear O D E : ( x \ ( T.\ ( ā11 ā T? A- T.\ y w X y W V -y x + x x - y 3 + y 2x2 -2z J (6.3) (6.4) and the in i t i a l region is the hypercube w i t h the d iagonal [ (1 .0 , -0 .05,0.9) , (1.2,0.05,1.1)] Because the mode l is nonl inear, we start by approx imat ing it w i t h an l inear differential inc lus ion. B y Tay lor 's theorem,V:rĀ®2/<g> z G [xi0,Xhi]Ā®{yi0,yhi]Ā®[zio: ZM}, the funct ion can be approx imated as: / x ^ ( ^ \ / x \ ( T. ā V J f y W y V 1 / X V X ā X y - y z ā ~z (6.5) where J is the Jacobianan Matrix and x = (xiā + Xhi)/2, y ā (yi0 + yhi)/2, z = (zio + Zhi)/2. A f te r comput ing the Jacob ian M a t r i x , the equivalent converts to: y ( x \ y \ * J ' 1 - 3x2 - 1 1 1 - 3y2 0 K Ax 0 - 2 J 71 0 \ / x \ ( 2x2 \ y \ z J + 2y2 -2z2 (6.6) Now, let us f ind the error bound : I Ix 3 - 3 x 2 x + 2 x 3 ^ A / = |/ - / I = (6.7) \y3 -3y2y + 2y3\ VI 2 ( x - x ) 2 \J Let A x = x ā x, and we rewri te x 3 - 3 x 2 x + 2 x 3 as 3 x ( A x ) 2 + ( A x ) 3 . N o t i n g that A x ā¬ , ^f t ' " 1 ' 0 ] , i t follows that the m a x i m u m value of | x 3 - 3 x 2 x + 2 x 3 | occurs for A x = x f t ' ~ x ' Ā° s i gn (x ) . A t th is point | x 3 - 3 x 2 x + 2 x 3 | = A x ( A x + 3|x|). S imi la r reasoning appl ies to the other two terms of A / , and we conclude: m a x ( A / ) < A x 2 . ( A x + 3|x|) Ay2{Ay + 3\y\) \ V (6.8) J 2 Ax2 where A x = (xhi - x ; 0 ) / 2 , A y = [yhi - yi0)/2, Az = (zhi - zlo)/2. T h i s y ie lds the differential inc lus ion: / x \ / 1 - 3 x 2 - 1 y c 0 \ / x \ I 2x2 \ 1 1 - 3 y 2 0 4x 0 - 2 J y V 2 / + +cube / A x 2 ( A x + 3|x|) \ Ay2(Ay + 3\y\) V 2 A : (6.9) C O H O uses th is approx imated dynamic mode l to compute the reachable state space. F igures 6.2 and 6.3 show the project ions onto the x-y and x-z planes. Note that in bo th project ions, the region at the end is contained in the in i t ia l region. T h i s establishes the invar iance of the computed region. A l l t rajector ies f rom this invar iant set remain in the set for al l t ime. T h e example is new; thus it is ha rd to f ind other results w i t h which to compare our result. We t ry to compare our result w i t h the two d imensional V a n der P o l osci l lator example i n [GM99] . T h e models for var iables x and y are the same, therefore we can compare the project ion of our reachable space onto the x-y plane w i t h the result f rom Greenstreet99a. F r o m F igu re 6.2, it is easy to see 72 F igu re 6.3: T h e Resu l t of 3 V d P E x a m p l e : y-z p lane 73 that the trace converges much faster. In [GM99] , the reachable set is computed on ly for the upper half cycle; the lower half cycle is the same w i t h the above half by symmetry . However, the reachable set for the whole cycle is computed i n our example. A l t h o u g h it is not necessary to verify the safety propert ies, i t demonstrates that our a lgor i thm and implementat ion are much more robust that the or ig inal version of C O H O . T h e use of symmet ry in the or ig inal version was because the over-approx imat ion of the reachable space became excessive as the projectagon went around the subsequent corners. Fur thermore, the or ig inal version was not able to analyse the three d imensional osci l lator ( [GM99] presents the more s tandard two-dimensional version) because of numer ica l s tabi l i ty issues. We have solved a l l of these problems. 6.2 Experimental Result for L P solver and L P project T h e S ink and 3D V a n der P o l osci l lator examples require a large number of l inear program solut ions and po lygon project ions. For example, there are 87148 l inear programs solved to ta l ly in the Sink example and 2051893 l inear programs solved in the 3D V a n der P o l osci l lator example w i th more than 1000 t ime advance steps. T h e success of ver i fy ing these examples demonstrates the soundness of our l inear program solver and pro ject ion a lgor i thms, the efficiency of the pro jectagon representat ion of reachable states and the tightness of our over approx imat ion f rom the l inear program solver, the pro ject ion a lgor i thm and model ing l inear izat ion. F i r s t , we examine the error (the w id th of the interval for the cost) of the op t ima l value produced by the l inear program solver. A s seen in F igu re 6.4, most of results have relat ive error less than 1 0 - 1 4 , and almost a l l op t ima l results have error bounds less than 1 0 - 7 . W e observe that the relat ive error is approx imate ly bounded by the \Julp [unit of least precision), A double precis ion u lp is 2 ~ 5 3 wh ich is roughly 1 0 ~ 1 6 . ) . F igu re 6.4 shows the d is t r ibu t ion of relat ive error of the op t ima l results for bo th examples. 74 -16 -15 -14 -13 -12 -11 -10 -9 -8 -7 L log10(error) J -16 -15 -14 -13 -12 -11 -10 -9 -8 -7 L log (error) J Figure 6.4: E r r o r D i s t r i bu t i on of O p t i m a l Resul t of C O H O L inea r P r o g r a m Solver Nex t , we examine the pivot selection method that we presented i n C h a p -ter 4.2.1. In the S ink example, each l inear program is solved w i t h about 4.67 pivots, whi le the number is much larger in the 3 V d P example, it is about 8.3 p ivots per l inear program. T h i s is because the 3 V d P example is three d imens iona l and has more constraints for each l inear program; therefore, there are many more feasible vertices to v is i t before reaching the op t ima l point . We observe that most pivots have a unique more favorable basis, wh ich demonstrates the efficiency of our pivot selection rules. F igu re 6.5 shows the d is t r ibu t ion of the number of branches per pivot . In the S ink example, roughly 85 percent of pivots have a unique successor, and almost a l l the rema in ing 15 percent have on ly two branches. In the 3 V d P exam-ple, more t han 86 percent of pivots have a clearly unique b ranch ; about 12 percent 75 of pivots have two branches; less than 1 percent of pivots have three branches; about 0.5 percent of pivots have four branches; and it is rarely that a pivot has more than four branches. x 10 Sink * 194760810 12,96 2 3 4 Number of Branches per Pivot 3VdP 10 t 5 159090 93444 2 3 4 Number of Branches per Pivot 12,78 Figure 6.5: D is t r i bu t i on N u m b e r of Branches of a P i v o t T h e analysis of i l l -condi t ioned system exceptions shows the necessity of the hand l ing w i t h bad ly condi t ioned problems. There are 8692 i l l -condi t ioned l inear systems found in the S ink example, about 0.01 per l inear p rogram. T h e 3 V d P is more complex and therefore has more i l l -condi t ioned system except ions. There are about 0.75 mi l l i on exceptions total ly, wh ich is 0.36 per l inear p rogram. T h e exper imenta l da ta shows our est imate of the condi t ion number to be qui te effective. A l l except ions i n the S ink example and more than 99.7 percent of except ions i n the 3 V d P example are caught by the est imated condi t ion number. T h e method only 76 fails for less than 0.3 percent in the complex 3 V d P example. M o s t of the bad ly condi t ioned l inear systems have huge condi t ion number; more t han 81 percent have condi t ion numbers greater than 1 0 1 6 . T h e d is t r ibu t ion is shown in F igu re 6.6. Sink 6 7 8 ,9 10 11 12 13, 14 15 16 17 L log1Q(condition number) J F igu re 6.6: D i s t r i bu t i on of Cond i t i on N u m b e r for I l l -Cond i t ioned Prob lems A s descr ibed in Chap te r 4.5, we t r ied two methods for f ind ing an in i t ia l feasible basis. T h e quick a lgor i thm is performed first, and the "b ig M " method is used i f the quick a lgor i thm fails. However, the exper imenta l da ta shows the percentage of l inear programs for wh ich the quick a lgor i thm succeeds is quite low. There are on ly 2 percent of l inear programs of 3 V d P example that contain a s imple feasible basis; whi le it is a l i t t le more than 0.2 percent in the S ink example. It is reasonable that the 3 V d P example has more s imple feasible bases because has many more constraints than the Sink example; therefore, the probab i l i t y to f ind a 77 simple feasible basis is much higher. For higher d imens ional models, the method might show its advantages. However, for the examples considered here, there is no signif icant advantage to inc lud ing the "quick" method . 8000 7000 6000 0& - * ā number of linear programs per step, 0 number of exceptions per step 100 200 300 400 Advancing Step 500 600 F igu re 6 . 7 : N u m b e r of L inear Programs and Excep t ions per Step Typ ica l l y , the designer specifies a simple projectagon such as a hypercube for the in i t ia l reachable set. A s reachabi l i ty computa t ion progresses, C O H O produces more compl icated projectagons to approximate the curvature of the actua l t rajecto-ries. If the number of po lygon vertices grows at each step, then C O H O wou ld become imprac t ica l l y slow. T h e analysis of the number of l inear programs and except ions along advanc-ing steps shows that C O H O l imi ts the complex i ty of the projectagons quite success-fully. A s shown i n F igu re 6 . 7 , the number of l inear programs oscil lates around a 7 8 constant value, wh ich shows the computat ion complex i ty for each step remains rel-at ively constant. However, there are some steps that have a much larger number of l inear programs. T h i s is caused by the project ion a lgor i thm, as descr ibed i n C h a p -ter 5.2. Back t rack ing is performed when the a lgor i thm j umps to the next quadrant incorrect ly. In rare cases, many l inear programs must be solved to compute the cor-rect result. W e should improve the project ion a lgor i thm in the future. T h e number of except ions has s imi la r t rends, thus the number of except ions per l inear programs almost remains near ly constant. W i t h a work ing version of C O H O , we now have the oppor tun i t y to examine a large sample of the i l l -condi t ioned problems that arise in C O H O ' S analysis. T h i s should allow iis to test var ious conjectures as to their cause and may enable a more efficient so lut ion to these problems i n the future. In summary , the exper imenta l results have shown the efficiency and robust-ness of our C O H O implementat ion. 79 Chapter 7 Conclusion 7.1 What has been Accomplished We have implemented a new l inear program solver for use in the C O H O ver i f icat ion system and demonstrated that it is efficient, accurate and robust . T h e framework of th is solver was set-up by L a z a , who explored the specia l s t ructure of C O H O l inear programs, presented an efficient 0(n) l inear system solver and combined it w i t h the S imp lex a lgor i thm. We implemented and improved his l inear system solver and integrated it w i t h interval ar i thmet ic wh ich guarantees the over approx imat ion as required by ver i f icat ion. T h e a lgor i thm to handle bad ly condi t ioned problems finds a good approx imat ion of the op t ima l result and great ly improves the accuracy of the l inear program solver. T h e exper imenta l results demonstrate the accuracy, efficiency and robustness of our solver. T h e pro ject ion prob lem, comput ing the project ion of a projectagon onto a two-d imensional subspace, is solved by our project ion a lgor i thm using the C O H O l inear program solver. T h e a lgor i thm produces a projected po lygon that might be a slight over approx imat ion . We der ived an analy t ica l upper bound for th is error and showed f rom our exper iments that it is negligible i n pract ice. We re- implemented the port ions of C O H O that rely on l inear p rogramming to use our new algor i thms for so lv ing l inear programs and pro ject ing the feasible region 80 onto two-d imensional subspaces. T h e new version of C O H O has been appl ied to two pract ica l examples. T h e success of these ver i f icat ion demonstrates the correctness of the ideas of C O H O sys tem such as projectagon representat ion of reachable region, approx imat ion of the mode l , etc. T h e exper imenta l da ta f rom these examples further demonstrate the soundness of the l inear program solver and pro ject ion a lgor i thm; otherwise, the error f rom these funct ions.would have prevented C O H O f rom ver i fy ing these complex examples. In summary , th is research has solved the previously known problems of C O H O l inear programs and pro ject ion of a projectagon; fur thermore, the implementat ion has been shown to be efficient and accurate enough for pract ica l ver i f icat ion tasks. In the process, the code base for C O H O has been extensively revised and re-organized to faci l i tate further enhancements. We sketch some of these possibi l i t ies i n the next sect ion. 7.2 Suggestions for Further Research T h e s tudy of C O H O l inear programs and project ion of the projectagon is not yet complete. There are many problems to explore further. T h e C O H O ver i f icat ion too l has just begun to show its abi l i ty for ver i fy ing real c i rcui ts and other designs; th is has opened up many areas for further research. O u r imp lementa t ion of interval ar i thmet ic as descr ibed in Chap te r 3.2 over approx imates the error bound . T h e I E E E f loat ing point s tandard [Gol91] defines the relat ive error of f loat ing point computat ion to be less than ^ulp. However, C O H O overapproximates it by a ful l ulp. We could get a t ighter b o u n d and improve the speed of the software by using special ized round ing modes that are inc luded in the I E E E s tandard . In par t icu lar , the I E E E s tandard defines several round ing modes inc lud ing round-up and round-down. T h e upper bound and lower bound can be obta ined by these two rounding modes. However, the Java language does not suppor t it. T h u s , C or assembly language wou ld be needed to access these 81 features, and J N I (the " J a v a Nat ive Interface) could be used to integrate th is w i t h i n our Java packages. B y implement ing interval versions of some basic l inear algebra operat ions such as vector dot product , we could avoid the overhead of large numbers of calls between Java and C and realize much of the performance of a fu l ly opt imized implementa t ion . T h e op t ima l bases returned by C O H O l inear p rogram solver are on ly numer-ica l ly possible op t ima l ; as descr ibed in Chap te r 4.2.2, the bases re turned by the solver can inc lude infeasible, super-opt imal bases as wel l . T h e inc lus ion of ext ra bases is caused by the est imat ion of computat ion error f rom interval ar i thmet ic ; it could be avoided if there is no round off error, i.e., the number has arb i t ra ry pre-cis ion. A r b i t r a r y precis ion ra t iona l ( A P R ) ar i thmet ic is an al ternat ive to interval ar i thmet ic or f loat ing point computa t ion i n which a l l computat ions are exact. U s -ing A P R , we cou ld implement a l inear program solver that wou ld on ly re turn the op t ima l basis (or bases if the cost funct ion is or thogonal to a c r i t i ca l constraint) . U s i n g A P R , most of the issues descr ibed in Chap te r 4.2 wou ld be automat ica l ly addressed. W e can implement an arb i t ra ry precision ar i thmet ic package based on the Biglnteger class, or use a mathemat ica l l ib ra ry l ike G M P [The]. T h e disadvantage of A P R computat ion is that it is not d i rec t ly suppor ted by s tandard hardware. Thus , it w i l l run slower than the corresponding floating point or interval number methods. Fur thermore, A P R methods tend to use more memory, and for general computa t ion , this the size of an A P R tends to grow w i t h the length of the computa t ion . Fortunately, our implementat ion of S imp lex starts f rom the or ig ina l d a t a after each p ivot and should therefore avoid much of th is memory overhead. Fur thermore , a hyb r id a lgor i thm could be used to improve efficiency: the interval me thod wou ld be appl ied to f ind a smal l set of poss ib ly op t ima l bases i n the S imp lex a lgor i thm, the arb i t rary precision rat ional ar i thmet ic wou ld on ly be used at the last step to ident i fy the t ru ly op t ima l basis f rom th is set. T h e opt ima l i ty of a basis can be easi ly checked by comput ing its feasibi l i ty for p r ima l and dua l l inear 82 programs. Therefore, th is a lgor i thm should have efficienty that is comparable to that of the pure ly interval ar i thmet ic method. W i t h arb i t ra ry precis ion rat ional representat ion, computa t ion is error free even for bad ly i l l -condi t ioned problems. Therefore i l l -condi t ioned basis should no longer be a prob lem. For the hybr id method, once an i l l -condi t ioned prob lem is encountered, the solver could be restarted using A P R methods for the par t icu lar p rob lem. Therefore, the op t ima l basis is always avai lable, and it is efficient because the bad ly condi t ioned problems rarely happend. We don not have to wor ry about the op t ima l basis be ing t ru l y singular. T h i s wou ld mean that some constraint is exact ly imp l ied by some other set of constraints, and the singular basis describes a l ine. Because the feasible space is bounded by const ruct ion, there are other constraints that impose endpoints on the feasible por t ion of th is l ine. A basis that replaces the redundant constraint w i t h one of these endpoint constraints is bo th op t ima l and non-singular . For the error of our l inear program solver, we believe that the m a x i m u m error i n the cost is bounded by the diameter of the reachable regions t imes the square root of the machine precision. T h i s conjecture is suppor ted by our numer ica l da ta as shown in Chap te r 6.2, but we have not completed a formal proof. If A P R methods are used, the exact op t ima l basis can be found; therefore, the relat ive error in the op t ima l cost w i l l be due to the use of double precis ion numbers to describe the or ig ina l l inear program and to report the f inal result. Thus , the error w i l l be determined by the machine precis ion. Ano the r idea to improve efficiency is to accelerate the pivot select ion method. Cons ider how s implex does pivot selection: for each co lumn j , the marg ina l cost c] is computed by equat ion 4.6. C o m b i n i n g this w i t h equat ion 4.5, we ob ta in : cj = Cj-4-{A^A..d) (7.1) C o m p u t i n g tj by equat ion 4.5 takes 0(n) operat ions using our special l inear system solver. L ikewise, ca lcu lat ing takes 0(n) operat ions. T h u s th is step takes 0(n) 83 operat ions per co lumn considered. M a t r i x mu l t ip l i ca t ion is associative. Therefore, equat ion 7.1 can be rewr i t ten as: cj = Cj - ((^Ag1) ā¢ A : J (7.2) Now, we can compute ( c ^ A g 1 ) using our 0(n) l inear system solver. T h i s is the same for a l l co lumns considered so we only need to calculate it once. In other words, we first compute d = (<$Ag1) (7.3) and use it for each co lumn as Cj = Cj - d- A-,j (7.4) Here, we can take advantage of the sparsi ty of A . A n y given co lumn of A has either one or two non-zero elements. Thus , d ā¢ A:j can be calculated i n 0 ( 1 ) t ime. U s i n g this method , the computa t ion t ime for a pivot is reduced f rom 0(m(nā m)) to 0(n + m ) , where m is the number of constraints and n is the number of variables. In the current version of the project ion a lgor i thm, a l l l inear programs are treated as independent op t im iza t ion problems by the C O H O l inear p rogram solver. However, the sequence of l inear programs to f ind the vertices of the pro ject ion poly-gon has a specia l s t ructure. T h e op t ima l basis of such a l inear p rogram corresponds to a vertex of the polygon. T h e op t ima l bases of adjacent l inear programs corre-spond to the two endpoints of an edge. Therefore, these two op t ima l bases is on ly a single pivot away. T h e C O H O l inear program solver shou ld be able to exploi t th is special s t ructure and improve the efficiency of pro ject ion a lgor i thm greatly. For the C O H O system, obviously, th is is just a start . The re are many topics to s tudy in the future. T h e top pr ior i ty is to app ly it on more examples, especial ly non- l inear c i rcui ts for high-speed d ig i ta l appl icat ions. W e are also interested in examples f rom hyb r id systems and control as wel l as models of b io logical process. 84 Projectagons provide an attract ive way to represent h igh d imens iona l objects. Pro jectagons can represent non-convex objects, and operat ions on projectagons are efficient i n terms of bo th t ime and space. E x p l o r i n g the issues of projectagon geom-etry is another area for future work. For example, the intersect ion of two or more projectagons is obta ined by the intersect ion of their pro ject ion polygons. O n the other hand , projectagons are closed under neither complement nor un ion , and we do not know how to calculate the smallest overapprox imat ion of the un ion . Such operat ions wou ld be useful for implement ing frontier based methods for reachabi l i ty that could offer signif icant improvements of efficiency. L ikewise, we do not know of an efficient a lgor i thm for determin ing whether or not a projectagon is empty. It appears that C O H O is wel l -sui ted for a h igh ly para l le l imp lementa t ion . W i t h i n each t ime step, the operat ions for each edge of each face are independent; accordingly, they could be computed concurrent ly. T h e amount of d a t a output for each edge is smal l compared w i t h the amount of work done to advance it th rough a t ime step. Therefore, paral le l comput ing or cluster compu t i ng techniques should be appl icable for improv ing the speed of C O H O ver i f icat ion. T h e procedure for l inear iz ing the non-l inear O D E mode l as descr ibed i n Chap te r 2.2 is done manua l l y for each ind iv idua l examples. T h i s is too cumber-some for use by pract ic ing hardware designers or contro l engineers. T h e abi l i ty to approx imate the mode l automat ica l ly should be added in a future version. We are present ly work ing on this for c i rcui t ver i f icat ion appl icat ions. Once we have manua l l y constructed models for transistors, we believe that these models can be composed d i rect ly f rom the structure of the circui t to produce a mode l for the cir-cui t . T h u s , no further manua l l inear izat ion should be needed once we complete a generic mode l for t ransistors. 85 i Bibliography [ADM02] Eugene Asarin, Thao Dang, and Oded Maler. The d/dt tool for verification of hybrid systems. In Proceedings of the Fourteenth Conference on Computer Aided Verification, pages 365-370, Copenhagen, July 2002. Springer. [BA] Ell iot Joel Berk and C. Scott Ananian. JLex: A lexical analyzer generator for Java(TM). http://www.cs.princeton.edu/~appel/modern/java/JLex/. [Bla77] R . G . Bland. New finite pivoting rules for the simplex method. Math. Oper. Res., 2:103-107, 1977. [BMP99] Olivier Bournez, Oded Maler, and Ami r Pneuli . Orthogonal polyhedra: Repre-sentation and computation. In Proceedings of the Second International Work-shop on Hybrid Systems: Computation and Control, pages 46-60. Springer, 1999. L N C S 1569. [B079] J . L . Bentley and T. Ottmann. Algorithms for reporting and counting geometric intersections. IEEE Transactions on Computers, 28:643-647,1979. [BTOO] O. Botchkarev and S. Tripakis. Verification of hybrid systems with linear dif-ferential inclusions using ellipsoidal approximations. Hybrid Systems: Compu-tation and Control, pages 73-88, 2000. L N C S 1790. [Dan63] George B . Dantzig. Linear Programming and Extensions. Princeton University Press, 1963. [DM98] Thao Dang and Oded Maler. Reachability analysis v ia face lifting. In Thomas A . Henzinger and Shankar Sastry, editors, Proceedings of the First International Workshop on Hybrid Systems: Computation and Control, pages 96-109, Berke-ley, California, A p r i l 1998. 86 [DMn05] Marc Daumas, Guillaume Melquiond, and Cesar M u noz. Guaranteed proofs using interval arithmetic. In Proceedings of the 17th IEEE Symposium on Com-puter Arithmetic, pages 188-195, 2005. [Eme97] E . Allen Emerson. Model checking and the mu-calculus. Proceeding of the DIM ACS Symposium on Descriptive Complexity and Finite Model, 1997. [Pre05] Goran Frehse. PHAVer : Algorithmic verification of hybrid systems past HyTech. In Proceedings of the Fifth International Workshop on Hybrid Systems: Com-putation and Control, pages 258-273. Springer-Verlag, 2005. L N C S 3414. [GM98] Mark R. Greenstreet and Ian Mitchell . Integrating projections. In Thomas A . Henzinger and Shankar Sastry, editors, Proceedings of the First International Workshop on Hybrid Systems: Computation and Control, pages 159-174, Berke-ley, California, A p r i l 1998. [GM99] Mark R. Greenstreet and Ian Mitchell. Reachability analysis using polygonal projections. In Proceedings of the Second International Workshop on Hybrid Systems: Computation and Control, pages 103-116, Berg en Dal , The Nether-lands, March 1999. Springer. L N C S 1569. [Gol91] David Goldberg. What every computer scientist should know about floating-point arithmetic. ACM Computing Surveys, 23(l):5-48, March 1991. [GV96] Gene H . Golub and Charles F . Van Loan. Matrix Computations. Johns Hopkins, 3 r d edition, 1996. [Hen96] Thomas A . Henzinger. The theory of hybrid automata. In Proceeding of the Eleventh Annual IEEE Symposium on Logic in Computer Science (LICS'96), pages 278-292, 1996. [HHWT97] Thomas A . Henzinger, Pei-Hsin Ho, and Howard Wong-Toi. HyTech: A model checker for hybrid systems. Internation Journal of Software Tools for Technology Transfer, 1(1-2):110-122, December 1997.' [HPWT01] T . A . Henzinger, Joerg Preussig, and Howard Wong-Toi. Some lessons from the HyTech experience. In Proceedings of the 40th Annual Conference on Decision and Control, pages 2887-2892. I E E E Press, 2001. 87 [HS74] Morris W . Hirsch and Stephen Smale. Differential Equations, Dynamical Sys-tems, and Linear Algebra. Academic Press, San Diego, C A , 1974. [LazOl] Marius D. Laza. A robust linear program solver for projectahedra. Master's , thesis, Department of Computer Science, University of Bri t ish Columbia, Van-couver, B C , December 2001. [PS82] Christos H . Papadimitriou and Kenneth Steiglitz. Combinatorial Optimization: Algorithms and Complexity. Prentice Hal l , Englewood Cliffs, N J , 1982. [PS85] Franco P. Preparata and Michael I. Shamos. Computational Geometry: An Introduction. Texts and Monographs in Computer Science. Springer, 1985. [RLB01] J . Douglas Faires Richard L . Burden. Numerical Analysis. Seventh edition edition, 2001. [SHA] Frank Flannery Scott Hudson and C. Scott Ananian. C U P : L A L R parser gen-erator for Java. http://www.cs.princeton.edu/~appel/modern/java/CUP/. [SKOO] B. I . Silva and B . H . Krogh. Formal verification of hybrid system using Check-Mate: A case study. In American Control Conference, volume 3, pages 1679-1683, June 2000. ' [SK03] Olaf Stursberg and Bruce H . Krogh. Efficient representation and computa-tion of reachable sets for hybrid systems. In Proceedings of the Sixth Interna-tional Workshop on Hybrid Systems: Computation and Control, pages 482-497. Springer, 2003. [SR+00] B. I . Silva, K . Richeson, et al. Modeling and verifying hybrid dynamical sys-tems using checkmate. In Proceedings of the 4th International Conference on Automation of Mixed Processes (ADPM 2000), pages 323-328, 2000. [The] The G N U Project. G M P : Arithmetic without limitations. http://www.swox.com/gmp/. [VanOl] Robert J . Vanderbei. Linear Programming: Foundations and Extensions. Springer, second edition, 2001. [Win87] Wayne L . Winston. Operations Research: Applications and Algorithms. Indiana University, third edition, 1987. 88
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Coho : a verification tool for circuit verification...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Coho : a verification tool for circuit verification by reachability analysis Yan, Chao 2006
pdf
Notice for Google Chrome users:
If you are having trouble viewing or searching the PDF with Google Chrome, please download it here instead.
If you are having trouble viewing or searching the PDF with Google Chrome, please download it here instead.
Page Metadata
Item Metadata
Title | Coho : a verification tool for circuit verification by reachability analysis |
Alternate Title | Coho : a verification tool for circuit verification by reachability analysis |
Creator |
Yan, Chao |
Publisher | University of British Columbia |
Date Issued | 2006 |
Description | COHO is a verification tool for systems modeled by nonlinear ordinary differential equations (ODEs). Verification is performed using reachability analysis. The reachable space is represented by projectagons which are the polyhedron described by their projection onto two dimensional subspace. COHO approximates nonlinear ODEs with linear differential inclusions to compute the reachable state. We re-engineer the COHO system to provide a robust implementation with a clear interface and well-organized structure. We demonstrate the soundness and robustness of our approach by applying it to several examples, including a three dimensional, non-linear systems for which previous version of COHO failed due to numerical stability problems. The correctness of COHO strongly depends on the accuracy, efficiency and robustness of the linear program solver that is used throughout the analysis. Our COHO linear program solver is implemented by integrating the Simplex algorithm with interval arithmetic based on the framework set up by Laza [Laz01]. For highly ill-conditioned problems, an approximation that is very close to the optimal result is computed by our algorithm. This results in a very small error bound that allows COHO to successfully verify interesting examples. This thesis also presents an algorithm to solve the problem of projecting the feasible region of a linear program onto two dimensional subspaces. This algorithm uses the COHO linear program solver for efficiency and accuracy. We derive an analytical upper bound for the error and present experimental results to show that the errors are negligible in practice. |
Genre |
Thesis/Dissertation |
Type |
Text |
Language | eng |
Date Available | 2010-01-16 |
Provider | Vancouver : University of British Columbia Library |
Rights | For non-commercial purposes only, such as research, private study and education. Additional conditions apply, see Terms of Use https://open.library.ubc.ca/terms_of_use. |
DOI | 10.14288/1.0051509 |
URI | http://hdl.handle.net/2429/18181 |
Degree |
Master of Science - MSc |
Program |
Computer Science |
Affiliation |
Science, Faculty of Computer Science, Department of |
Degree Grantor | University of British Columbia |
GraduationDate | 2006-11 |
Campus |
UBCV |
Scholarly Level | Graduate |
AggregatedSourceRepository | DSpace |
Download
- Media
- 831-ubc_2006-0724.pdf [ 3.84MB ]
- Metadata
- JSON: 831-1.0051509.json
- JSON-LD: 831-1.0051509-ld.json
- RDF/XML (Pretty): 831-1.0051509-rdf.xml
- RDF/JSON: 831-1.0051509-rdf.json
- Turtle: 831-1.0051509-turtle.txt
- N-Triples: 831-1.0051509-rdf-ntriples.txt
- Original Record: 831-1.0051509-source.json
- Full Text
- 831-1.0051509-fulltext.txt
- Citation
- 831-1.0051509.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
data-media="{[{embed.selectedMedia}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
https://iiif.library.ubc.ca/presentation/dsp.831.1-0051509/manifest