UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

Formal verification of a digital PLL Wei, Jijie 2014

Your browser doesn't seem to have a PDF viewer, please download the PDF to view this item.

Item Metadata

Download

Media
24-ubc_2014_september_wei_jijie.pdf [ 3.88MB ]
Metadata
JSON: 24-1.0165974.json
JSON-LD: 24-1.0165974-ld.json
RDF/XML (Pretty): 24-1.0165974-rdf.xml
RDF/JSON: 24-1.0165974-rdf.json
Turtle: 24-1.0165974-turtle.txt
N-Triples: 24-1.0165974-rdf-ntriples.txt
Original Record: 24-1.0165974-source.json
Full Text
24-1.0165974-fulltext.txt
Citation
24-1.0165974.ris

Full Text

Formal Verification of A Digital PLLbyJijie WeiB. Computer Science, Harbin Normal University, 2009M. Computer Science, Harbin Engineering University, 2012A THESIS SUBMITTED IN PARTIAL FULFILLMENTOF THE REQUIREMENTS FOR THE DEGREE OFMaster of ScienceinTHE FACULTY OF GRADUATE AND POSTDOCTORALSTUDIES(Computer Science)The University Of British Columbia(Vancouver)August 2014c© Jijie Wei, 2014AbstractCommon AMS circuit are composed from blocks that can be modeled accuratelyusing linear differential inclusions to enable verification of important properties us-ing reachability analysis. This dissertation presents a formal verification of DigitalPhase Locked Loop (PLL) using reachability techniques.PLLs are ubiquitous in analog mixed signal (AMS) designs and are widely usedin modern communication equipment, clock generation for CPUs in computers,clock-acquisition in high-speed links etc. The most important property of a PLLis convergence, which means starting from any possible initial conditions, the PLLwill eventually lock to the desired equilibrium. We model the digital PLL as aset of Ordinary Differential equation (ODEs), and discretize the weakly non-linearODEs to linear differential inclusions. The transformation not only provides us anover approximation of the verification problem but also provides the basis for asound proof.We present the verification of a digital PLL using real tools, SPACEEX andCOHO-REACH. In particular, we show how each component of the digital PLLcan be modelled as a hybrid automaton. Due to the large number of transitionscaused by the model, the whole proof is established by several lemmas. Interest-ing problems such as a timing glitch in the Phase Frequency Detector (PFD) arediscussed and triggering conditions are formally proved in the dissertation. Globalconvergence is demonstrated by both tools.Based on the digital PLL circuit and the challenges that arose during our ver-ification, the error bounds, limitations, implementation differences and usabilityof the two leading tools are carefully evaluated. SPACEEX provides a graphicaluser interface that makes it easy to get started with simple examples but requiresiiextensive user interaction for larger problems. The interface to COHO-REACH is aMATLAB API. While this lacks the packaged-tool feel of SPACEEX, it provides aflexible way to break a large verification problem into smaller lemmas and allowsthe proof to be “re-executed” simply by re-executing the MATLAB script.iiiPrefaceChapter 3. Figure 3.5 is provided by Yu Ge with permission. Portions of text areused and modified from [65] of which I am one of the authors. Prof. Greenstreetbuilt the behavior model for the digital phase locked loop (PLL) from [20]. I wrotethe MATLAB code for simulation and proved the PFD glitch safe bound.Chapter 4. A version of the material has been published as J. Wei, Y. Peng, G. Yuand M. Greenstreet. (2013)[65]. Two sets of parameters (Set 2. and Set 3.) fromTable 4.1 are provided by the author of the design in [2, 20]. I performed all theexperiments and established the global convergence using SPACEEX.Chapter 5. Prof. Greenstreet derived the linear model (Equation B.9) for lemma 4.I performed the experiments and established the global convergence using COHO-REACH.ivTable of ContentsAbstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iiPreface . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ivTable of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vList of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viiiList of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ixAcknowledgments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xi1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11.1 AMS Circuits and Verification . . . . . . . . . . . . . . . . . . . 11.2 Hybrid System Verification . . . . . . . . . . . . . . . . . . . . . 41.3 Phase Locked Loop Verification . . . . . . . . . . . . . . . . . . 51.4 Contributions . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71.5 Overview . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72 Related Work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 92.1 Formal Verification of AMS Circuits . . . . . . . . . . . . . . . . 92.1.1 Equivalence Checking . . . . . . . . . . . . . . . . . . . 102.1.2 Model Checking . . . . . . . . . . . . . . . . . . . . . . 112.1.3 Deductive Methods . . . . . . . . . . . . . . . . . . . . . 122.2 Formal Verification of Hybrid Systems . . . . . . . . . . . . . . 132.2.1 Hybrid Automata . . . . . . . . . . . . . . . . . . . . . . 13v2.2.2 Reachability Analysis Tools . . . . . . . . . . . . . . . . 142.3 Formal Verification of Phase Locked Loops . . . . . . . . . . . . 192.4 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 203 Digital Phase Locked Loop . . . . . . . . . . . . . . . . . . . . . . . 223.1 Phase Locked Loop . . . . . . . . . . . . . . . . . . . . . . . . . 223.2 A Digital PLL . . . . . . . . . . . . . . . . . . . . . . . . . . . . 253.2.1 A Digital PLL . . . . . . . . . . . . . . . . . . . . . . . 263.2.2 Behavior Model . . . . . . . . . . . . . . . . . . . . . . 333.2.3 PFD Glitch . . . . . . . . . . . . . . . . . . . . . . . . . 363.3 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 394 Verification of Global Convergence using SPACEEX . . . . . . . . . 404.1 Ideal SPACEEX Model and Explosion . . . . . . . . . . . . . . . 404.1.1 Ideal SPACEEX Model . . . . . . . . . . . . . . . . . . . 414.1.2 Explosions . . . . . . . . . . . . . . . . . . . . . . . . . 434.2 Global Convergence Verification . . . . . . . . . . . . . . . . . . 434.2.1 Model 1lo . . . . . . . . . . . . . . . . . . . . . . . . . . 464.2.2 Model 2lo . . . . . . . . . . . . . . . . . . . . . . . . . . 484.3 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 505 Verification of Global Convergence using COHO-REACH . . . . . . 515.1 Introduction and Motivation . . . . . . . . . . . . . . . . . . . . 515.2 Global Convergence Verification . . . . . . . . . . . . . . . . . . 555.2.1 Initial Trajectories . . . . . . . . . . . . . . . . . . . . . 565.2.2 Leaving the Wall . . . . . . . . . . . . . . . . . . . . . . 585.2.3 Convergence of fDCO to fref . . . . . . . . . . . . . . . . 605.2.4 Convergence of c to ccenter . . . . . . . . . . . . . . . . . 615.3 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 656 Conclusion and Future Work . . . . . . . . . . . . . . . . . . . . . . 666.1 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 666.1.1 Global Converge Verification by SPACEEX . . . . . . . . 676.1.2 Global Converge Verification by COHO-REACH . . . . . . 67vi6.1.3 Pros and Cons of SPACEEX and COHO-REACH . . . . . . 686.2 Future Work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72A PFD Glitch Proof . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 80A.1 Proof . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 80B Model Derivation for fDCO and Z3 Proof . . . . . . . . . . . . . . . . 87B.1 Z3 Validation of the Derived Model . . . . . . . . . . . . . . . . 90C Model Derivation for Spaceex Model 2 . . . . . . . . . . . . . . . . . 96D Coho Verification Code . . . . . . . . . . . . . . . . . . . . . . . . . 99viiList of TablesTable 2.1 Comparison of Reachability Analysis Tools . . . . . . . . . . 15Table 3.1 Four Sets of Parameters . . . . . . . . . . . . . . . . . . . . . 34Table 3.2 Four Sets of v Bounds for PFD Glitch . . . . . . . . . . . . . 38Table 4.1 Four Sets of Parameters . . . . . . . . . . . . . . . . . . . . . 45Table A.1 Four Sets of Parameters . . . . . . . . . . . . . . . . . . . . . 84viiiList of FiguresFigure 1.1 A Digital Radio Receiver . . . . . . . . . . . . . . . . . . . . 2Figure 1.2 Typical AMS Circuit and Corresponding Hybrid System . . . 4Figure 1.3 A General PLL Block Diagram . . . . . . . . . . . . . . . . 6Figure 3.1 A Simple PLL . . . . . . . . . . . . . . . . . . . . . . . . . 23Figure 3.2 A Simple Low Pass Filter . . . . . . . . . . . . . . . . . . . . 25Figure 3.3 A Simple Charge-Pump PLL . . . . . . . . . . . . . . . . . . 26Figure 3.4 The Digital PLL Verified . . . . . . . . . . . . . . . . . . . . 27Figure 3.5 Ring-Oscillator Response . . . . . . . . . . . . . . . . . . . 28Figure 3.6 Digitally-Controlled Oscillator . . . . . . . . . . . . . . . . . 29Figure 3.7 The PFD . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31Figure 3.8 A Example of Harmonic Lock for XOR PD . . . . . . . . . . 32Figure 3.9 Trajectory . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35Figure 3.10 Trajectories . . . . . . . . . . . . . . . . . . . . . . . . . . . 35Figure 3.11 PFD Glitch . . . . . . . . . . . . . . . . . . . . . . . . . . . 37Figure 4.1 Ideal Hybrid Automata for Each Component of the Digital PLL 42Figure 4.2 SPACEEX Simulation of Limit Cycle Oscillation for Variable c 44Figure 4.3 Hybrid Automaton for the PFD . . . . . . . . . . . . . . . . 46Figure 4.4 Convergence of v . . . . . . . . . . . . . . . . . . . . . . . . 49Figure 5.1 MATLAB Simulation of One Step of the Limit Cycle Oscil-lation, (a) Expected, (b) Support Function with Four Directions 52ixFigure 5.2 SPACEEX Simulation of One Step of the Limit Cycle Oscilla-tion, (a) uses Support Function with 512 Directions, (b) usesDefault Settings in SPACEEX . . . . . . . . . . . . . . . . . . 53Figure 5.3 COHO-REACH Modes and Switching Surfaces . . . . . . . . 57Figure 6.1 SPACEEX Simulation, (a) uses Default Settings in SPACEEX,(b) uses Support Function with 512 Directions . . . . . . . . 69Figure A.1 Change of c in Region One, Four and Seven . . . . . . . . . . 85Figure B.1 Variable Shifting . . . . . . . . . . . . . . . . . . . . . . . . 88xAcknowledgmentsFirst, I would like to express my sincere gratitude to my supervisor, Dr. MarkGreenstreet, for his steady support of my Master’s studies and research, for hispatience, motivation, enthusiasm, and wisdom. He introduced to me the researcharea of formal verification, and helped me throughout the research and writing ofthis thesis. In addition, Dr. Mark Greenstreet is a great mentor and friend, whoalso gave me advice on how to do better research and better programming. I couldnot have imagined a better supervisor for my Master’s study.Besides, I would like to thank my seconder reader, Prof. Ian Mitchell, for hisinsightful comments and valuable questions. I am grateful to Prof. Shahriar Mirab-basi for his precious help with Phased Locked Loop.Many thanks to Chao Yan for his support on COHO-REACH and discussions onthe verification of the Phase Locked Loop. Also, I would like to thank my fellowlab mates in ISD group: Yan Peng, Shabab Hossain, Brad Bingham, Sam Bayless,Celina Val, Mike Enescu for the stimulating discussion and the fun we had forthe past two years. And I greatly appreciate the help from my friend Ge Yu onSpectre R© simulation.Last but not least, I would like to thank my family, my parents and sister, forgiving me unconditional support and encouragement.xiChapter 1IntroductionAnalog-mixed-signal (AMS) design has emerged as the dominant paradigm forimplementing analog functionality in deep-submicron integrated circuit technolo-gies. These designs consist of a small number of analog blocks with a mixture ofdigital and analog control loops to ensure correct operation. Phase Locked Loops(PLLs) are one of these widely used feedback loop systems. Such systems can bemodelled as hybrid automata; then reachability techniques can be used to validateproperties such as global convergence, the focus of this dissertation. In this chap-ter, Section 1.1 provides a brief overview of analog mixed signal (AMS) circuitsand the verification challenges they pose. AMS circuits can be naturally modelledas hybrid systems, and Section 1.2 presents the main idea behind hybrid systemsand associated verification methods. Section 1.3 combines these to look at the spe-cific AMS that we study in this dissertation: Phase Locked Loop. Section 1.4 andSection 1.5 summarize the contribution of this thesis.1.1 AMS Circuits and VerificationElectronic circuits can be divided into two broad categories: analog and digital.Most integrated circuits today are primarily digital. Digital designs benefit fromsophisticated automatic design tools, and digital circuits can operate correctly evenwith the far-from-ideal properties of very small transistors. However, the real worldis continuous, if a digital electronic device is to interact with real world, it will al-1LPFLPFADCADCPre-ampAntennaFigure 1.1: A Digital Radio Receiverways need an analog interface. For instance, a digital radio receiver (see Figure 1.1)has an analog front end typically consisting of a pre-amplifier, a mixer to convertthe signal to baseband or an intermediate frequency, possibly more amplifiers, ademodulator, and an analog to digital converter. The design of analog circuits ismuch less automated than the design process for digital. The main tools are circuitsimulators that are still very much like the original SPICE [56]. Furthermore, ana-log circuits are much more sensitive to the properties of individual transistors anddo not operate well with the low power supply voltages required by modern digitalcircuits. For these reasons, analog circuits are being replaced by AMS designs. AnAMS circuit consists of both digital and analog components where digital blocksare used to replace the functions that have traditionally been implemented withanalog and to provide control and adaptation functions that enable the remaininganalog circuits to operate in spite of the challenges presented by circuit character-istics in state-of-the-art fabrication processes.Examples of mixed-signal integrated circuits include data converters using delta-sigma modulation [14], analog to digital converter [60], digital radio chips [16] anddigital phase locked loops [21].Most ICs produced currently use CMOS technology. Emphasis on high levelsof system integration creates a high demand for ICs that are predominantly digitalbut include many AMS subsystems. In addition, the explosion in wireless commu-nications, such as cell phones, tablets, and the emerging “internet-of-things”, mo-tivates the widespread use of AMS designs. Significant numbers of mixed-signal2chips are also used in set-top boxes, multimedia and consumer products such ascamcorders, and industrial devices such as temperature sensors, medical devices,and industrial controllers. Industry analysts have estimated that 80% of all IC de-signs today are mixed signals [11].As the use of AMS in modern system on chip (SoC) designs increases, theverification tasks become significantly more complex and time-consuming. Today,formal methods such as equivalence checking and model checking are a standardpart of digital design with more advanced approaches used for especially criticalfunctions. In contrast, analog circuit verification remains a difficult research prob-lem. Analog signals are a natural part of all continuous applications, includingaudio, radio, images and sensor for temperature, motion, etc. Analog circuits pro-vide the interface between these real-world inputs and outputs and the digital worldof processors and memory. Although analog circuits comprise a small fraction ofthe total circuitry on most chips, they account for a large fraction of the designerrors [66]. The combination of the discrete dynamics and continuous dynamicsin AMS design further complicates the verification problem, which requires a highlevel of expertise and deep understanding of the circuit’s behavior.The traditional approach of validating AMS circuit is simulation using SPICE,Spectre, or similar tools. While it works for DC and operating point analysis,simulation is inadequate for verifying critical global properties such as convergenceto the desired operating point from all initial conditions. Even one simulation froma single initial point could take weeks, months, or longer, and getting reasonableconverge of the range of possible starting condition is impractical. In addition,simulation falls short in validating more interesting properties of the design suchas temporal requirements. For instance, the lock time of Phase Locked Loop.The success of formal methods for digital designs verification has motivatedexploring applications of formal approaches to other domains. AMS designs can benaturally modeled as hybrid systems: state machines with a continuous dynamicsgiven by ordinary differential equations (ODEs) at each state, and transitions whoseguard conditions depend on the values of continuous variables. This motivates ourinterest in hybrid systems which are introduced in the next section.3ActuatorSensorDigital ComputerPlantDACADCDigital CircuitAnalogCircuitA typical hybrid system A typical AMS circuitFigure 1.2: Typical AMS Circuit and Corresponding Hybrid System1.2 Hybrid System VerificationHybrid systems [47] combine both discrete and continues dynamics. An typicalexample of hybrid systems is modelling a system where a digital computer is con-trolling a continuous “plant” such as a chemical manufacturing process or an auto-mated vehicle. AMS circuits can be modelled as hybrid systems where the digitalcircuit is receiving quantized inputs from the analog blocks, and outputs of the dig-ital circuit set modes and serve as inputs to the analog block. Figure 1.2 showsa typical AMS circuit and a typical hybrid system. Extensive research has beendone in modelling and verification of the hybrid system. This thesis focuses onapplications of reachability analysis techniques.Reachability analysis for hybrid systems is a challenging and active researcharea. There are a wide range of verification techniques and tools for reachabilityanalysis. State-of-the-art computational tools for hybrid system are generally oftwo kinds: those that use symbolic computation to calculate reachable sets whilelimiting the continuous dynamics to simple abstractions such as rectangular inclu-sions; and tools that use optimization techniques to approximate the set of reach-able states by polyhedra [19], ellipsoids[52], level set [57] etc.Tools such as HYTECH [44–46] and DREAL [34] use symbolic computationmethods. On the other hand, D/DT [8, 22], CHECKMATE [17, 61], and PHAVER[30, 31] use optimization techniques. The tools used in this dissertation are SPACEEX[33] and COHO-REACH [37, 67], which are also based on optimization. Althoughthey both use polyhedra to represent the reachable states, SPACEEX uses piecewise4linear models that are given by a fixed partitioning of the state space called a “hy-bridization”. This causes SPACEEX to incur mode transitions as an artifact of thispartitioning. COHO-REACH, in contrast, performs on-the-fly linearization whichcan be done for the entire reachable region if the non-linear terms in the model aresmall, or on a face-by-face basis for more highly non-linear models. This allowsCOHO-REACH to have smaller over-approximations and to avoid the overhead ofextra mode transitions. SPACEEX and COHO-REACH are described in more detailin Section 2.2.2.1.3 Phase Locked Loop VerificationPhase-locked loops (PLLs) are ubiquitous in analog and mixed signal (AMS) de-signs. Phase locked loops are used as frequency multipliers to generate the multi-GHz clocks for CPUs from a crystal oscillator reference (typically around 50−100MHz). By making the multiplication factor programmable, CPUs can switchbetween different clock frequencies to obtain peak performance on-demand andsave power when there is less work to do. PLLs are also used to coordinate thetiming on high-speed interfaces such as DDR memory, Ethernet interfaces, andUSB. Other applications include modulators and demodulators for wireless com-munication as well as generating the carrier frequency for radio transmitters andreceivers.A phase locked loop is a closed-loop frequency control system. A typical PLL(see Figure 1.3) circuit consists of phase detector (PD), low pass filter, voltage-controlled oscillator (VCO) and divider. The phase detector compares the referencesignal and output signal of the frequency divider, generating an output that is ameasure of their phase difference. The phase error is used to adjust the VCO sothat the phase of the VCO phase locks to the phase of the input signal. The dividerproduce an output signal whose frequency is 1/N that of the VCO. Thus, whenthe PLL adjusts the VCO frequency so that the phase and frequency of the divideroutput match that of the reference, then the VCO frequency will be N times that ofthe reference, and their phases will be aligned.A failed PLL can block further test of an entire chip or major subsystems;thus, there is a high value in verifying the correctness of PLL designs. Functional5 	!"	Figure 1.3: A General PLL Block Diagramverification of analog blocks such as PLLs can be divided into two parts: globalconvergence to an intended operating point and small signal analysis at the operat-ing point [49]. The key insight here is that nearly all analog blocks are intended tohave some kind of linear response when at or near their intended operating point.Existing analysis techniques such as periodic AC analysis (PAC) [50] are availablein standard commercial CAD tools such as Spectr R© from Cadence. These tech-niques allow designers to characterize key performance properties of analog blockssuch as the jitter, power-supply sensitivity, and tracking bandwidth of a PLL. A de-signer can have high confidence in the correct functioning of their block assumingthat it reaches its intended operating point.To show that an analog block will reach its intended operating point from anyinitial state is the global convergence problem. Here, the non-linearities of thecircuit must be taken into account. There are basically two approaches for handlingthe non-linearities; one is to discretize the non-linear part, so that the model ispiece-wise linear [3, 4]; the other is to employ symbolic techniques [27, 48, 55,64]. Simulation based methods are impractical both because of the impossibilityof covering all possible initial points in a continuous space and because individualsimulations are very compute intensive. Verification of global convergence wouldbe of value both to find bugs in existing designs and to enable more aggressive PLLdesign. Section 2.3 surveys work on verifying global convergence of PLLs.61.4 ContributionsCommon AMS circuits are composed from blocks that can be modeled accuratelyusing linear differential inclusions to enable verification of important properties us-ing reachability analysis.The research described in this dissertation presents formalverification of a digital phase locked loop using reachability techniques. Specifi-cally, there are three main contributions:1. The behavior model of the digital PLL is represented as a set of ODEs. Ourfirst contribution is the transformation of the weakly non-linear ODEs anddiscrete time recurrences to linear differential inclusions. This ensures thatthe reachability tools compute an over approximation of the reachable spaceand provides the basis for a sound verification.2. We present the verification of the digital PLL using real tools. In particu-lar, we establish global convergence from any initial state the PLL eventu-ally reaches a state of phase and frequency lock using SPACEEX [33] andCOHO-REACH [68]. We show how each component of the digital PLL canbe modeled as a hybrid automaton. Our models account for non-linearitiesof the components, quantization, and other non-idealities. We demonstratehow global convergence can be shown by reachability analysis.3. Based on the digital PLL circuit and the problems arose during our verifica-tion work, we evaluate the advantages, differences, both algorithmic and inimplementation of these two leading tools for hybrid systems when appliedto AMS circuit verification.The basic components of oscillators, dividers, phase comparators, and inte-grators are common to all PLL designs, We expect this work to be re-usable forverifying other PLL structures. We also believe that the linear differential inclu-sion approach and the verification methods used here are applicable to other AMScircuit verifications.1.5 OverviewThis dissertation is organized as follows:7• Chapter 2 describes prior research on AMS verification, PLL verification andhybrid systems. The mathematical representation of the hybrid system is for-mally presented in this section. In addition, Chapter 2 explores reachabilityanalysis techniques and well developed tools.• Chapter 3 presents the digital PLL that we verified. It describes the behav-ioral model of the circuit as well as its dynamics. The property to be verifiedand derivation of the non-linear differential equations are formally presentedin this chapter.• Chapter 4 describes the modelling of the digital PLL as a hybrid automatausing SPACEEX. Appendix B shows how we obtain a linear differentialinclusion (LDI) from the non-linear ODE of Chapter 3. Section 4.2 presentsour verification of the global convergence of the digital PLL, that is, fromany valid initial state, the PLL will converge to the desired region. Theadvantages and limitations of SPACEEX are discussed.• Chapter 5 presents the modelling of the digital PLL as hybrid automata us-ing COHO-REACH. Different from SPACEEX verification, we use COHO-REACH as a lemma prover and establish the global convergence by severallemmas.• Chapter 6 concludes the dissertation. SPACEEX and COHO-REACH are com-pared based on the verification of the digital PLL. Ideas for future works arepresented.8Chapter 2Related WorkThis section provides a survey of prior work on AMS circuit verification. Sec-tion 2.1 gives an overview of formal methods in AMS circuits verification, whichincludes equivalence checking, model checking and deductive methods. Notingthat AMS circuits fall naturally in the category of hybrid system, Section 2.2presents prior works on hybrid system verification. Hybrid systems are often mod-elled using hybrid automata which are presented in Section 2.2.1. In addition,a wide range of reachability analysis tools are discussed and compared in Sec-tion 2.2.2. Finally, Section 2.3 describes previous work on Phase Locked Loopverification.2.1 Formal Verification of AMS CircuitsThe verification of AMS designs is concerned with the assurance of correct func-tionality, for instance, convergence and metastability verification. Other propertiesof interest include whether an AMS design is robust with respect to different typesof uncertainties like the affect of temperature, jitter, etc. Traditionally, SPICE sim-ulation has been used in critical operating points for verification where the test isoften done manually in an informal fashion. This barely qualifies as verificationas the search of the state space is not complete. Like pure digital circuits, formalmethods for AMS circuits can be grouped into three classes: equivalence checking,model checking and reachability analysis, and deductive methods. To cover the9complete state space, most of these methods use over approximation techniques.This section reviews the research works in the field of AMS design formal verifi-cation.2.1.1 Equivalence CheckingEquivalence checking is the problem of showing that two system models are equiv-alent with respect to some notion of functional similarity in terms of input-outputbehavior. The two models to be compared may be either at the same levels of ab-straction or at different levels, for example, equivalence of behavioral models andimplemented transistor net list or equivalence of two behavioral models. It allowsthe user to find, analyze and eliminate all the errors introduced during transfer fromone abstraction to another.Equivalence checking for digital circuits has achieved great success early onand has been widely used. Generally, the abstraction of the circuits are formallypresented as boolean formulas and there are two basic techniques used in equiva-lence checking:1. Binary decision diagrams (BDDs) [15]: A specialized data structure de-signed to support reasoning about boolean functions.2. Conjunctive Normal Form Satisfiability: SAT [53, 62] solvers returns anassignment to the variables of propositional formula that satisfies the formulaif such an assignment exists.Equivalence checking for analog circuits is more challenging. As the statespace of analog circuits is continuous, transformation or discretization is needed.In [13], the authors applied equivalence checking between the implementationand specification of the circuit, which are represented by transfer functions andare transformed into the discrete Z-domain. By doing this, the transfer functionsare digitized and encoded into finite state machine (FSM) representations such asBDDs. As the transfer function generation for non-linear systems is impractical,this methodology is only suitable for linear systems.In [42], the author proposed an equivalence checking approach between twonon-linear analog circuits. The system is described by sampling the state and com-10puting differences in the vector fields with the aid of non-linear transformations(Jordan form) on the sampled space.Although, there are many kinds of transformations, still it is non-trivial work,and the error caused by discretization needs to be carefully included to make theverification sound. In addition, most of the equivalence checking work on AMScircuit separates digital components and analog components, for example in [59]the coupling between the analog and digital parts are ignored.2.1.2 Model CheckingModel checking was first proposed in [18, 58] as a paradigm for computer-aidedverification of finite-state programs and has gained widespread acceptance. Givena model of a system, model checking tools exhaustively and automatically checkwhether this model meets a given specification. To this end, both the model andthe specification are formulated in some precise mathematical language, such as afinite-state machine derived from a hardware description language program for themodel and a temporal logic formula for the specification. Often, the specificationis for safety properties that must hold in all states; such properties can be written asassertions, i.e. propositional logic formulas. The approach is especially appropriatefor the design and verification of digital circuits and distributed protocols.The first work extending model checking techniques to AMS systems was byKurshan and McMillan in [51]. The authors partitioned the continuous space intoa fixed grid of hyper rectangles. The reachability relations between these cubeswere computed using numerical techniques. This idea was then generalized andimplemented in the tool AMCHECK [41], which also took into account delaytimes during discretization of the state space.Representing the reachable space by enumerating hyper rectangles results ina computation whose time and space requirements grows exponentially with thedimensionality of the continuous space. To overcome the expensive computation,Greenstreet and Mitchell [37] use projectahedra (use discretization and projectiontechniques.) to reduce dimension. While the approach is less precise, over ap-proximations are systematically employed to make sure that the verification is stillsound. This idea was further developed and implemented in tool COHO-REACH by11Yan [67]. COHO-REACH has been used to verify AMS circuits including a high-speed toggle [69], a Mo¨bius Ring-Oscillator [71], an asynchronous, handshakingarbiter [70] and more recently a digital phase locked loop [65].Variant approaches of polyhedral based analysis and related tools such as HYTECH[44–46], D/DT [8, 22], CHECKMATE [17, 61], PHAVER [30, 31], SPACEEX [33]etc. are summarized in Section 2.2. Gupta et al. [40] utilize CheckMate to ver-ify analog circuits such as a tunnel diode oscillator and a delta-sigma modulator.Dang et al. [23] use D/DT to verify a biquad low-pass filter. Frehse et al. [32]used PHAVER to verify analog oscillator circuits. Wei and Peng et al. [65] usedSPACEEX to verify the global convergence of digital PLL. the whole validation isdivided into verification of several lemmas as SPACEEX is not able to handle thewhole system at a time mostly due to its over approximations.Most of the existing model checking tools can handle small or linear systems.When the system become more complex or includes some non-linear parts, the run-ning time maybe unaffordable. Thus, dividing the verification into the validation ofa set of smaller lemmas is a practical choice which requires a deep understandingof the dynamics of the circuit to be verified, in other words, white-box verification.2.1.3 Deductive MethodsAutomatic theorem proving (ATP or automated deduction) proves mathematicaltheorems based on inference rules by computer programs or with the aid of a userdevised deductive proof. In practice, ATP methods often require large efforts byverification experts, and the theorems that are proven tend to be of a prosaic natureto show hardware or software correct rather than being mathematically profound.1There have been several attempts to use ATP to verify properties of AMS circuits.In [35] Gosh et al. used the PVS theorem prover verify the equivalence ofpiece-wise linear synthesized analog circuits and their user given behavioral spec-ification. The synthesized circuits model is extracted from the sized componentnet list produced by the synthesis tool, in terms of characteristic behavior of thecomponents and various voltage and current laws. The verification was appliedfor DC and small signal analysis. It was a good starting point of a new method-1An example of ATP for proofs that are mathematically interesting: http://www.cse.chalmers.se/research/group/logic/TypesSS05/Extra/wiedijk 2.pdf12ology for AMS circuit design, the non-linearity was handled by piece wise linearapproximation, yet no formal error analysis was included.In [1], the authors proposed a new method to prove the properties of AMSdesigns. In this approach, the AMS circuit is represented by a system of recurrenceequations (called “SREs in paper). These equations model the dynamics of theanalog components sampled at clock cycles of the digital part. Properties of theAMS circuit are proven using Induction. The symbolic simulation and the proofstrategy are implemented inside the computer algebra system Mathmatica, whichprovides special functions to simplify and prove algebraic relations. The processis iterative and is able to generate counterexamples once the proof failed. Thoughthe methodology can be automatized and integrated into the design flow, additionalwork is needed in order to generate the SRE model from circuit descriptions. Theauthors illustrated their method by analyzing a third-order, delta-sigma, analog-to-digital converter.In [65] Wei and Peng use Z3, a SMT solver, to verify the global convergenceof a simplified model of a digital PLL, the validation of the convergence propertywas transformed into the validation of Lyapunov function.2.2 Formal Verification of Hybrid SystemsHybrid system theory lies at the intersection of the areas of engineering controltheory and formal verification. The continuous and discrete dynamics in hybridsystem not only coexist, but also interact and change in response to discrete events.The continuous dynamic are described by a set of differential equations or inclu-sions. Significant research has been done in modelling and verification of hybridsystem.2.2.1 Hybrid AutomataA wide range of models for hybrid systems has been proposed. The most com-monly used models for hybrid system include hybrid automata, transition systems(or switched systems) and hybrid petri nets. The common features of all of thesemodelling paradigms is the interaction of different dynamics.A hybrid automata [43] is a finite state machine with a finite set of continuous13variables whose values are described by a set of ordinary differential equations. Itconsists of modes and the transitions between those modes. Each mode is associ-ated with a continuous dynamics. For example, these dynamics can be described asdifferential equations such as x˙ = f (x), where x ∈Rn evolves continuously in time.The mode transitions allow the system to change dynamics when the continuousstate satisfies a specific guard transition. These transitions are treated as occuringinstantaneously. Furthermore, each mode has an invariants, a condition that musthold for any trajectory in that mode, transition must occur before the invariant isviolated.There are many variants of hybrid automata: timed automata (TA) [6] which in-clude a “clock” measure with derivative 1, and linear hybrid automata (LHA) [43]with linear dynamics.In other words, each mode defines a polyhedral constraint (not necessarily con-vex) for x˙. The constraints for x˙ are independent of x. The invariant, µ and guardsδ are defined using polyhedra in a similar way.2.2.2 Reachability Analysis ToolsHybrid system verification is an active research area. Problems such as modeswitching in autopilot, and air traffic management (ATM) can be modeled as hy-brid automata. Safety problems are critical in these problems, especially whenhumans are involved. Other interesting verification properties include stability,convergence and robustness. Taking AMS circuit for instance, the verification ofcritical circuit parts is able to find bugs before the design is fabricated as a chip.Verification of these circuits not only saves billions of money but also shortens thetime to market. Thus, the application of formal verification methods for hybridsystems to AMS circuits is of great, pratical importantance.Reachability analysis takes the models of hybrid systems and explores the statespace completely by solving both the continuous and discrete dynamics. It is apromising techniques and widely used in hybrid system verification. Based ondifferent kinds of models, there are a wide range of state-of-the-art reachabilitytools. Table 2.1 summarizes and compares the feature of a set of reachability toolswith the remainder of this section giving a more detailed description of each tool.14Table 2.1: Comparison of Reachability Analysis ToolsTool Model Representation Linear NonlinearHYTECH LHA symbolic,convexpolytopesconstant ODIs2 hybridizationD/DT LDHA orthogonalpolyhedronlinear ODIs hybridization,face liftingCHECKMATE PIHA, DTTS flow pipe linear ODEs flow pipeCOHO-REACH NHA projectagon linear ODIs hybridizationPHAVER linear HIOA convexpolyhedronconstant ODIs hybridizationSPACEEX LDHA supportfunctionsconstant ODIs piecewiseaffineDREAL NHA Interval linear ODEs nonlinearODEsHYTECH , developed by Henzinger et al. [44–46], is one of the earliest toolsfor hybrid systems modeled by linear hybrid automata (LHA). It uses a symbolicmodel checking approach and is able to analyze a set of concurrent automata andperform parametric analysis. HYTECH computes the condition under which a lin-ear hybrid system satisfies a temporal requirement. HYTECH is able to generatesa diagnostic error trace if the verification fails. It should be noted that HYTECHcannot handle hybrid systems with linear differential equations because the con-straints on the flows (i.e. x˙) depend only on the mode and are independent of thevalues of the continuous variables (i.e. x) within the mode.D/DT [8, 22] is a reachability analysis tool developed originally by Thao Dang.It computes the reachable sets for models with linear continuous dynamics, e.g. x˙ =Ax+b. The reachable sets are represented by orthogonal polyhedra, i.e. polyhedrawhich are finite unions of full-dimensional hyper-rectangles. The computation isbased on solving linear Ordinary Differential inclusions (ODIs) using the maxi-mum principle method [63] and nonlinear ODEs using face-lifting, a method torepresent non-convex polyhedra. Later extensions include support of hybridiza-tion [9, 10] and CEGAR [5, 7]. D/DT has been used for safety verification. How-ever, D/DT has a problem that the number of hypercubes grows exponentially withthe dimension of the problem.2ODIs: Ordinary Differential Inclusions15CHECKMATE [17, 61] was developed at Caregie Mellon University. It is aMatlab-based tool for simulation and verification of threshold-event-driven hybridsystems (TEDHS). CHECKMATE represent sets as convex polyhedra and computesthe reachable states by propagating polyhedra under linear and affine dynamics,which, similar as D/DT over approximates nonlinear dynamics along each surfaceof the polyhedra based on Pontryagin Maximum Principle [63].COHO-REACH was developed by Yan, Greenstreet and Mitchell [37, 67] forcircuit verification. It represents reachable sets using two-dimensional projectionsof higher dimensional non-convex polyhedra and evolves these “projectagons” un-der affine over approximations of nonlinear dynamics using linear programming.More details of COHO-REACH are presented later in this section.PHAVER [30, 31] is a tool for verifying safety properties of hybrid system.It was targeted to overcome the limitation of HYTECH , which suffers from over-flow due to the use of exact arithmetic with a data structure of limited digits. Tothis end, PHAVER use Parma Polyhedra Library (PPL) [12] and the GNU MultiplePrecision Arithmetic Library (GMP) [36] to remove these numerical restrictions. Itrepresent sets as polyhedra. PHAVER ’s model is based on linear automata, whichare defined by linear predicates and piecewise linear constant, convex, bounds onderivatives. Affine dynamics are handled by on-the-fly over approximation. Alsofrom the Verimag lab, SPACEEX [33] was later developed to facilitate the imple-mentation of algorithms related to reachability and safety verification. More detailsof SPACEEX are presented in later in this section.More recently, Gao and Clarke et al. [34] developed an SMT solver for nonlin-ear formulas over the reals call DREAL . Similarly to HYTECH , it uses a symbolicmodel checking approach. DREAL uses Interval Constraint Propagation (ICP) asits main framework to find solutions of real constraints. While the framework hasbeen proposed earlier on in [28, 29], DREAL is based on δ complete decision pro-cedure, where δ is the user specified error bound: it returns either unsat or δ -saton input formulas.16COHO-REACHCOHO was original proposed by Greenstreet and Mitchell [37] and developed byYan and Greenstreet [67]. It was originally optimized for the highly nonlineardynamics of VLSI circuits modeled at the transistor level. Later, an extensionof COHO , COHO-REACH, was developed that added support for mode switch-ing along with optimizations to improve the performance and reduce the over-approximation error for weakly non-linear models.COHO-REACH was developed and implemented for a simple class of moderate-dimensional hybrid systems with nonlinear ODE dynamics. It represents reach-able sets using projectagons. Projectagons provide a compact representation ofhigh-dimensional geometric object by projecting them onto two-dimensional sub-spaces. These two-dimensional projections are computed using linear program-ming. COHO-REACH solves nonlinear ODEs by approximating ODEs as lineardifferential inclusions. This construction involves a matrix exponentiation thatCoho-Reach bounds by using the Taylor expansion of the exponential along witha bound for the tail of the expansion. COHO-REACH constructs systems of lin-ear inequalities that bound the solutions to the differential inclusions. Coho-Reachmakes extensive use of linear programming and checks its solutions using intervalarithmetic and the “complementary slackness” condition. If these indicate that acomputation may be too inexact, Coho-Reach repeats the computation using ar-bitrary precision rational arithmetic to obtain exact solutions. COHO-REACH isoptimized with a range of techniques, including a hybrid application of double-precision, interval, and arbitrary precision rational arithmetic and various approx-imation algorithms with provable error bounds to reduce error as well as improveperformance. COHO-REACH has been used to verify AMS circuits like High-Speed Toggle [69], Mobius Ring-Oscillator [71], Arbiter [70].SPACEEXSPACEEX [33] is a platform for verifying hybrid system. The goal is to verify thata given mathematical model of a hybrid system satisfies the desired safety prop-erties. It combines both polyhedra and support function representations of contin-uous sets to compute an over approximation of the reachable states. SPACEEX’s17analysis core is not a single tool but a development platform on which several dif-ferent verification algorithms are implemented. Each algorithm may implement itsown set representation and apply to its own class of models. Such a bundle of analgorithm, set representation, and model is referred to as a scenario. Currently, twoscenarios are implemented:• The PHAVER scenario implements the basic algorithm from the tool PHAVER[30, 31]. It applies to Linear Hybrid Automata, which are hybrid systemswith piecewise constant bounds on the derivatives.• The LGG Support Function scenario implements a variant of the Le Guernic-Girard (LGG) algorithm [39]. It applies to hybrid systems with piecewiseaffine dynamics with nondeterministic inputs.There are two main improvements compared to PHAVER . First, the approxi-mation model is based on a first-order approximation of the interpolation error andthe accuracy is further improved by combining the forward and backward approx-imation. Second, SPACEEX employs variable time steps to guarantee a given errorbound.COHO-REACH and SPACEEX are the two verification tools used in this disser-tation. There are three main differences between the approach taken by SPACEEXand that as implemented in COHO-REACH. First, SPACEEX uses hybridization[10, 24] to model nonlinear systems. Hybridization, as implemented in SPACEEX,statically partitions the state space into small regions and bounds the nonlinear dy-namics by simple differential inclusions in each region. COHO-REACH performson-the-fly linearization to construct inclusions as they are needed. This producestighter inclusions and avoids switching events that are merely artifacts of the hy-bridization process and not actual switching behaviours of the underlying system.Second, SPACEEX represents reachable polyhedra using a set of support vectorsthat is fixed throughout the analysis. The SpaceEx user-interface lets the userchoose several pre-defined forms for the the support vectors. The most accurateof these is to use N vectors at angles of 0, 360/N, 2 ∗ (360/N) 3 ∗ (360/N) de-grees in each of the (d choose 2) planes formed by pairs of coordinate axes fora d-dimensional system. This amounts to projecting the reachable region onto18planes, much as has been done in Coho [38]. SpaceEx also allows the user to pro-vide a set of support function vectors. COHO-REACH determines support vectorsby adaptively by computing the time advanced normals of the polyhedrons faces.Third, SPACEEX uses a simple norm-based method to handle inclusions and finiteapproximations of matrix exponentials. COHO-REACH uses an error bound withdirection, and adapts the number of terms in its approximations to achieve a desiredtolerance.2.3 Formal Verification of Phase Locked LoopsPhase Locked Loops (PLLs) are a common AMS circuit (Section 1.1). PLLs arewidely employed in radio, telecommunications, computers and other electronicapplications. The function of a phase-locked loop (PLL) is to adjust the PLL’soscillator so that it tracks the frequency and phase of a reference signal. A criticalproperty of a PLL is lock at the right frequency.There have been several previously published reports of Phase locked loopverification using formal methods. The earliest verification that we know of was by[26]. Dhingra’s design uses a fixed frequency oscillator that is intended to operateat N times the frequency of a reference. The PLL adaptively chooses edges of itsinternal oscillator to approximate the edges of the lower frequency reference. Thetime resolution is limited by the period of the internal oscillator. While this maybe useful for low frequency applications such as audio frequency modems, we arenot aware of any such PLLs in use for the more standard PLL applications of clockgeneration, clock data recovery, and wireless communication. Dhingra verified thetracking behaviour of his design using the HOL theorem prover.More recently, Dong et al. [27, 64] proposed using property checking for AMSverification, including PLLs. They used “symbolic recurrence equations” as aproperty specification language and show how this can be used to automaticallyconstruct a monitor to check simulation runs to see if a PLL locks in the requiredtime for that run. This does not address the problems of long-simulation times toshow that a PLL locks or the incompleteness of simulation based approaches toshow convergence.Shortly after the work by Dong et al., Jesser and Hedrich [48] described a19model-checking result for an analog PLL with an XOR-gate phase detector. Theyperformed symbolic model checking using multi-terminal binary decision diagrams(MTBDDs) to represent both the discrete and analog parts. They state that thefour-dimensional analog state space is partitioned into 211 hyper boxes, and thatnext-box relations are determined by random simulation. It is not clear how theyguarantee the complete coverage with this approach.Later, Althoff et al. [3, 4] presented the verification of a charge-pump PLLusing an approach that they refer to as “continuization”. They use a purely lin-ear model for the components of their PLL, and their focus is on the switchingactivities of the phase-frequency detector, in particular, uncertainties in switchingdelays. Their approach verifies the PLL for ranges of component parameters.More recently, Lin et al. [54] independently developed an approach for verify-ing a digital PLL using SMT techniques. To the best of our knowledge, they are thefirst to claim formal verification of a digital PLL. Similar to the approach presentedin this paper, they consider a purely linear, analog model and then reason about thediscrepancies between this idealized model and a digital implementation. They usethe KRR SMT solver to verify bounds on this discrepancy. They verify bounds onthe lock time of a digitally intensive PLL assuming that most of the digital vari-ables are initialized to fixed values, and that only the oscillator phase is unknown.Our work shows initialization for a different PLL design over the complete statespace.2.4 SummaryThis chapter surveyed related research work on AMS circuit verification. Equiv-alence checking, model checking and deductive methods were described. Equiv-alence checking is used to compare the equivalence of two models either at thesame level of abstraction or at different levels of abstraction. BDD and SAT are thetwo most commonly used techniques for discrete circuit equivalence checking. Toapply equivalence checking to analog circuits, transformations are needed. Modelchecking checks the correctness of a property given a formally presented model.To this end, both the model and specification are formulated in a precise mathe-matical language. Automatic theorem proving for AMS circuit verification is still20a premature approach. Related tools includes PVS, Z3 etc.AMS circuits include both discrete and continuous parts, it is natural to modelthem as a hybrid system. Section 2.2 presents hybrid automata, one commonlyused model for hybrid systems. It also describes reachability analysis, a promisingverification method for state exploration in hybrid system. In addition, a widerange of reachability analysis tools, for example HYTECH, D/DT, COHO-REACHetc. are discussed and compared.A phased locked loop is a typical AMS circuit that is ubiquitous in SoC (System-on-Chip) designs and are of great importance because the failure of a PLL wouldcause huge impact. Phase locked loops are widely used in analog and mixed signaldesigns. Recently, PLL design has shifted from traditional, analog, charge-pumpbased designs to various “digital” architectures. Several consequences of devicescaling to smaller feature sizes have motivated this transition including: greaterdevice-to-device parameter variation impairs designs that depend on matched cir-cuits; lower power supply voltages removes the “voltage headroom” needed forhigh-quality, on chip current sources; and the scaling of passive components suchas inductors and capacitors lags that for transistors. Therefore, the coexistence ofboth analog and digital parts in a general PLL increases the difficulty of the veri-fication work. PLL verification is a critical topic for both researchers in academiaand engineers in industry. Section 2.3 reviewed the related research work on PLLverification. In the next chapter, we will present our digital PLL model and derivethe ordinary differential inclusions used for our global convergence verification.21Chapter 3Digital Phase Locked LoopPhase Locked Loops are versatile circuits that are ubiquitous in AMS designs.Section 3.1 presents the general structure of a phase locked loop. Section 3.2describes the digital PLL that we verified in this dissertation. More specifically,Section 3.2.1 explains the detailed components of the digital PLL. Section 3.2.2presents the behavioral model and analyzes the dynamics of the digital PLL. Forthe behavior model of the digital PLL, hybrid automata provide a natural model formost components of the digital PLL. However, the phase-frequency detector (PFD)has behaviors that are not naturally modelled using hybrid automata. Section 3.2.3describes these issues and presents my solution. The detailed mathematical proofof this solution presented in Appendix A.3.1 Phase Locked LoopThe function of a phase-locked loop (PLL) is to adjust the PLL’s oscillator so that ittracks the frequency and phase of a reference signal. This section sketches the mainprinciples of analog PLL designs as these motivate digital PLL designs, and mostof the function blocks of a digital PLL correspond to blocks of a traditional, analogPLL. Figure 3.1 shows a simple PLL. The PLL sets the control voltage, vctrl of theVCO according to a combination of the phase difference between the VCO and thereference and the integral of this difference. This matches the VCO’s frequency tothat of the reference and aligns their phases. Simply, if the PLL’s oscillator lags22LPFVCOPhase comparator∆!!!!"# !!!"# !!!"# !Figure 3.1: A Simple PLLthe reference, then vctrl increases; and the VCO frequency increases so that theVCO will catch up with the reference. Conversely, if the PLL’s oscillator is aheadof the reference, then vctrl will decrease causing the PLL’s oscillator frequency todecrease.In more detail, the voltage-controlled oscillator (VCO) can be understood asa voltage-to-frequency converter. Phase is the integral of frequency; so we canexpress the phase of the VCO output as θvco = (∫vctrldt) mod [−pi,+pi). Phase ismodular, and we write θ mod [−pi,+pi) to indicate the value in [−pi,+pi) that iscongruent with θ modulo 2pi radians1. The phase comparator generates an outputvoltage that is proportional to the phase difference of its inputs: θ∆ = (θvco−θre f )mod [−pi,+pi). The reference is assumed to have a constant frequency, ωre f ; thusθre f = ωre f t, where t is time. The low pass filter, LPF, implements the integraland proportional correction terms with vctrl = a0θ∆+ a1∫θ∆dt. Combining theseequations yields:θvco = a0∫θ∆dt +a1∫ ∫θ∆dtdt (3.1)We differentiate twice to obtain:d2dtθvco = a0ddtθ∆+a1θ∆ (3.2)In the simple case where a0 = 0 and a1 < 0, the PLL becomes a simple har-monic oscillator. The frequency of the PLL oscillates with mean value of ωre f butnever converges. If both a0 and a1 are negative, then there is a unique solutionwhere the PLL’s oscillator converges to the frequency and phase of the reference.Note that if all of the components are linear, then simple algebraic techniques suf-fice to show (or refute) global convergence.1If the phase detector has internal state, then the interval for ‘mod’ can be wider23The phase comparator from Figure 3.1 can be made from an exclusive-OR(XOR) logic gate. We interpret a logically true output of the XOR gate as +1and a logically false output as −1. The low-pass filter output includes the integralof the average value of the output. To converge, this average must be zero. Thismeans that the phase difference between the reference and VCO is 90 degrees.For stability, this average must be positive when the VCO lags the reference andnegative when the VCO leads. Thus, the PLL locks with the VCO lagging thereference by 90 degrees. The PLL converges when the difference of phase angle is90 degrees.A simple first order Low Pass Filter (LPF) consists of a resistor and a capac-itor. As shown in Figure 3.2, the input signal Vin is applied to the series combi-nation (both the resistor and the capacitor together) while the output signal Vout istaken across the capacitor only. Let Z1 be the impedance of the resistor, Z2 be theimpedance of the capacitor, ω be the radians of the input signal. We know thatVoutVin=Z2Z1 +Z2=1jωC1jωC +R||vout ||||vin||=1√1+(ωCR)2(3.3)Therefore, if ω  1RC and the magnitude of the gain is approximately unity, theoutput of equals the input. If ω  1RC then the gain goes to zero. Thus f =ω2pi =12piRC is called the cutoff frequency. Designers often want low loop bandwidth forstability and low jitter. This means that the capacitor must be large which can takeup an impractical amount of silicon area. Digital PLLs overcome these problemsby replacing the low pass filter with discrete-time digital filter.Another common design is a charge-pump PLL. A charge-pump PLL can bebuilt by simply connecting all the above components using a charge pump betweenthe phase-detector and the loop-filter as an integrator. Figure 3.3 shows this simplePLL. One source of the charge-pump is connected to the positive supply rail whilethe other is connected to the negative supply rail. The sources are separated by twoswitches S1 and S2. The output of the phase detector controls the gating signal.24RCVinVoutZ1Z2Figure 3.2: A Simple Low Pass FilterWhen the PD detects that fref is ahead of fDCO, S1 is on and S2 is off which causescurrent to flow out of the pump and into the LPF. On the other hand, when freflags behind fDCO, current flows out of LPF and into the pump. Charge pump PLLshave been popular for several reasons including that they lock with a 0-degreephase offset between the VCO and the reference and the loop gain and bandwidthcan be adjusted by adjusting the pump current IP.Various applications modify these basic designs. For example, a divider canbe included between the VCO and the phase comparator to obtain frequency mul-tiplication. A PLL can be used as a phase or frequency modulator by adding amodulation signal to the control voltage. Conversely, it can be used as a demodu-lator by outputting the control voltage.For real designs, the components are not perfectly linear. These componentsmay very closely match their linear idealizations when the PLL has locked to thereference, but significant non-linearities may occur when out of lock. Furthermore,the analog components of the simple PLL are difficult to implement in advancedfabrication technologies. For example, large capacitors are needed to implementthe integrator part of the low pass filter. For these and other reasons, designers aremaking more and more extensive use of digital and mixed signal designs for PLLs.3.2 A Digital PLLThe digital PLL verified in this dissertation is described in Section 3.2.1. Accord-ing to the dynamics of this digital PLL, the behavior model is presented as a set ofdifferential equations in Section 3.2.2. The digital PLL consists of a PFD which25VCOPDCharge-pump!!! !!!! !+!! −!! !! !! !!"# !!"# LPF!! !: charge-pump currentFigure 3.3: A Simple Charge-Pump PLLcompares the arrival times of edges from the (frequency divided) DCO and the ref-erence to generate the up or down signal; a BBPFD connected with a c integratorand a v integrator, which together control the digital VCO; a LPF, a ∆Σ modulatorand a digital Voltage-controlled Oscillator. While most of these components canbe naturally modelled by hybrid automata, the PFD has a potential clock-glitchissue that we were unable to represent using simple hybrid automata models. Theclock glitch in the PFD happens when one of the input signal goes up during areset, and the built-in delay would cause the up signal to be ignored. The problemis important for it could cause the DCO controls to head away from convergence.Section 3.2.3 discusses these issues.3.2.1 A Digital PLLFigure 3.4 shows the digital PLL architecture from [20]. The three control paths ofthe PLL (linear phase, bang-bang frequency, and coarse frequency) work togetherto set the frequency, fDCO, of the digitally controlled oscillator (DCO) to N timesthe reference frequency, fref , and to align the phase of the DCO and the reference(i.e., rising edges of the reference clock should coincide with rising edges of theDCO).26vΣFrefΣFrefDACCdecap−+CoarseControlFrequencyBang−BangFrequencyControlLinearPhaseControlFref FDCOFDCO0:34:7∆ΣBBPFD0:230:1415:23LPF0:7DCOPFD+−dnup∆θcCenter   code ÷NFigure 3.4: The Digital PLL VerifiedThe digitally-controlled oscillator in [20] is a three-stage ring-oscillator withthree control inputs: v, c, and ∆θ . The upper sub figure of Figure 3.6 shows onestage of the DCO. The applied input voltage determines the instantaneous oscil-lation frequency. In addition, each stage is shunted with a set of capacitors thatregulate the DCO frequency, the larger the capacitance connected to the stage, thelonger the period and the lower the frequency. In the design we analyzed, four ca-pacitors translated a digital value into the corresponding capacitance. The remain-ing capacitor is used by the ∆Σ modulator to interpolate between discrete settingsof the capacitance.The v input sets the operating voltage of the DCO. First-order transistor mod-els suggest that fDCO should be roughly proportional to v. Figure 3.5 (a) showsthe results of Spectre R© simulations (provided by Ge Yu [72]) of a ring-oscillatorin a 65nm CMOS process with a 1.2V nominal Vdd . For an operating voltage vwith 0.5 ≤ v ≤ 1.2, a linear fit provides a good approximation of the DCO fre-quency. Likewise, the c input controls switches that add capacitors to increase theload for each stage of the ring oscillator. As seen in Figure 3.5 (b), the oscillatorperiod (inverse frequency) is very accurately modeled as a linear function of loadcapacitance.Note that the DCO frequency, fDCO is proportional to a linear function of v and270.4 0.5 0.6 0.7 0.8 0.9 1 1.1 1.200.511.522.533.5operating voltage, v (volts)DCO frequency, f DCO (GHz)  00.511.522.533.5simulation datalinear fit(a) frequency vs. operating voltage1 2 3 4 5 6 7 8 9 1000.511.522.533.5load capacitance, c (pF)DCO frequency, f DCO (GHz)  00.511.522.533.5simulation datalinear fit for 1/fDCO(b) frequency vs. load capacitanceFigure 3.5: Ring-Oscillator Response28updnupdnaVbC7C6C5C4∆Σ One Stage of the DCOThree-Stage DCOVupdnab∆Σ Vupdnab∆Σ Vupdnab∆Σ c2c4c8c!!:!! !!:!! !!:!!Figure 3.6: Digitally-Controlled Oscillator29inversely proportional to a linear function of c. We can write this asfDCO(c,v) =1+αv1+βc f0 (3.4)where f0 is the frequency of the DCO (extrapolated) to where v and c are at zero,and α and β are the sensitivities to v and c respectively. Most importantly, whenmodeling the response to both v and c, either their ranges must be quite small, orthe model is inherently non-linear. To show global convergence, the non-linearityof the DCO must be included in the model. For simplicity we “normalize” v, c andfDCO in the remainder of this dissertation and simplify Equation 3.4 to :fDCO(c,v) = vc (3.5)The ∆θ input controls a linear phase path. Each stage of the ring-oscillatorconsists of an inverter in parallel with two tri-state inverters. One of the tri-stateinverters is enabled when up is asserted, and the other is disabled when dn isasserted. This causes the oscillator to run faster when up is asserted without dn,and slower when dn is asserted without up. This produces, a phase shift advance(resp. delay) that is proportional to the length of time that the oscillator is in its fast(resp. slow) mode. When both up and dn are asserted at the same time, the DCOfrequency is the same as when neither is asserted (to within the design tolerancesof the inverters and transmission gates).The frequency divider, labeled ÷N in Figure 3.4 is a modulo-N counter. Itsoutput has a period N times that of the DCO.The two phase-frequency detectors (PFD and BBPFD) in the diagram deter-mine if the DCO should increase or decrease its frequency. Figure 3.7 shows atypical circuit for the linear phase-frequency detector, labeled PFD in Figure 3.4,and its state transition diagram. Each state is labeled with the value of the up anddn (down) outputs of the PFD – for example, in state 01 the up signal is low andthe dn signal is high. The PFD is reset to state 00 after it has received rising edgesfrom both the reference clock, ref, and the frequency divided DCO, DCO/N. If thenext event is a rising edge from the divided DCO, this indicates that the DCO isahead of the reference in phase or running at a higher frequency. Accordingly,30DCO/N001001DCO/N11refDCO/NrefState DiagramD QRDQRdnupCircuitrefFigure 3.7: The PFDthe PFD enters state 01 (asserting dn) indicating that the DCO should decrease infrequency (equivalently, drop back in phase). Likewise, if the first event after re-setting is a rising edge of the reference, then up is asserted. Note that the durationof asserting dn without up (or vice-versa) indicates the phase difference betweenthe reference and DCO. This is the principle behind the linear-phase control pathof the PLL. If the DCO frequency is much higher than the reference, then DCO↑event will always precede the ref↑ event, so the PFD will signal “down” on eachcycle. Likewise, if the DCO frequency is much lower than the reference, the PFDwill signal “up” on each cycle.Phase-frequency detectors offer several advantages over phase-detectors forPLL designs. PDs are sensitive to input duty cycle and will lock with a phaseerror if the reference or DCO/N duty cycle is not 50%. PD based designs are alsosusceptible to harmonic lock. Harmonic lock occurs when the PLL locks with theDCO frequency (divided by N) is a multiple or sub multiple of the reference. Forinstance, when the reference frequency fDCO is twice of fref , and the phase differsin 90 degree. Figure 3.8 shows this case. The last signal (Φvco′ in the figure) whosefrequency is twice of fref , outputs a phase difference that is the same as fDCO (Φvcoin the figure). The PLL can be locked with either Φvco or Φvco′ . The first, Φvco is thedesired lock. The second, Φvco′ , is an example of harmonic lock —the frequencyof Φvco′ is twice the frequency of the reference.In contrast with an XOR PD, a PFD resets its internal state when both signalsare high and is not sensitive to the input duty cycle. Furthermore, PFD based31!!"# !!"#!!!"#!!!"#′!Figure 3.8: A Example of Harmonic Lock for XOR PDdesigns do not suffer from harmonic lock because the PFD responds to the “first”edge after reset of each input.The BBPFD is a so-called “bang-bang” phase frequency detector. It simply de-tects whether the linear PFD asserts up before dn or dn before up and then outputs+1 or −1 respectively. The output of the BBPFD indicates the direction that theDCO should move, but does not contain the magnitude information of the linearPFD.The accumulators implement integrators. The integrator in the bang-bang pathdirectly controls the c input of the DCO. The Delta-Sigma modulator provides av-eraging: if the lower four bits of the accumulator for c encode a value of k, then theoutput of the Delta-Sigma modulator will be asserted for k, roughly evenly spaced,out of 16 cycles of the DCO. The integrator in the coarse frequency control pathis designed to have low-bandwidth to ensure low jitter (cycle-to-cycle variation ofthe DCO period). The output of the integrator is converted to a voltage with thedigital-to-analog converter (DAC). To suppress the coupling of power-supply noiseinto the DCO, an additional low-pass filter and linear regulator are included in theloop.To understand the operation of the digital PLL when 1N fdco 6= fref , the differ-ence will be noted by the PFD and the BBPFD. The BBPFD will output a value todrive c in the direction to correct for the frequency difference. For small frequencydifferences, this bang-bang loop provides fast tracking. For larger differences, theaccumulator for c will saturate at its minimum or maximum value. Any of thesesituations lead to c being different from the ccenter. The PLL is designed so thatc provides fast-tracking in response to disturbances. By adjusting v to return c to32ccenter, the PLL maintains good tracking range. This drives the coarse frequencyloop to change v in the direction needed to bring fDCO/N to fref and to return c tothe center code value. Because the DCO has only a fixed set of values, for mostvalues of fref there is no choice of c and v that causes fDCO/N to be exactly fref .This leads to limit-cycle oscillations of the DCO control settings. The linear phasecontrol path reduces the magnitude of these oscillations. The linear phase path alsoprovides stability to the loop. Many of the design choices in [20] are made so thatthe contribution of the digital control to the jitter of the DCO output is less than theunavoidable jitter due to thermal and shot noise.3.2.2 Behavior ModelTo use verification tools for hybrid systems, we approximate the discrete dynamicsfor the DCO control capacitance, c, and voltage, v using ordinary differential equa-tions (ODE). Section 4.2.1 and Appendix B show how these differential equationscan be transformed into linear differential inclusions that include error terms thataccount for approximating the discrete dynamics with continuous ones, and forapproximating a non-linear system as a piece-wise linear one. For simplicity, weomit the Delta-Sigma modulator, low-pass filter, and voltage regulator from thismodel, and assume a divider ratio of 1. The ODE that we use is:satc = ((c = cmin)∧ (θδ < 0))∨((c = cmax)∧ (θδ > 0))satv = ((v = vmin)∧ (c > ccenter))∨((v = vmax)∧ (c < ccenter))Resetθ∆ = ((θ∆ ≤−2pi)∧ (θ∆ ≥ 2pi)c˙ = −sign(θ∆)g1 fref , if ¬satc= 0, if satcv˙ = g2 fref (c− ccenter), if ¬satv= 0, if satvθ˙∆ = fDCO(c,v)− fref −Kθ∆ if ¬Resetθ∆reset θ∆ to θ∆−2pisign(θ∆), if Resetθ∆(3.6)We describe each of these equations below.33c˙ is determined by the sign of the phase difference. When the phase differenceis positive, that is fDCO is ahead of fref , c will increase (g1 is negative and gives thesize of a “step” in c for a single cycle of fref ) until the phase difference changesits sign. Table 3.1 shows four sets of possible parameters. The sign of the phasedifference corresponds to the output of the BBPFD as is shown in Figure 3.7. TheBBPFD drives c to a value where fDCO ≈ fref . Then, the coarse voltage path setsv to a value so that c is near its centre value with fDCO nearly equal to fref . If theinitial value of fDCO(c,v) is much different than fref , then the accumulator for cwill saturate at cmin or cmax.Table 3.1: Four Sets of Parametersparameters Set 1. Set 2. Set 3. Set 4.c [0.9, 1.1] [0.9, 1.1] [0.9, 1.1] [0.9, 1.1]v [1.0, 2.5] [0.2, 3.8] [0.1, 1.5] [1.5, 2.5]θ 5.0 5.0 5.0 5.0fref 2.0 3.0 1.0 2.0ccenter 1.0 1.0 1.0 1.0g1 -0.01 -0.003125 -0.003125 -0.01g2 g1/5 g1/5 g1/5 g1/5K 0.1 0.8 0.8 0.7v is the output of the LPF according to Figure 3.4. For simplicity, we omit theLPF and voltage regulator in this model; thus, v is the output of the v-accumulator.v˙ changes with the difference of c and ccenter, which is the centre value for c. Fig-ure 3.9 shows two trajectories of the digital PLL model starting from two points inv× c space. The blue one shows when fDCO is initially much lower than fref . Thevalue of c decreases and saturates at the minimum point, which drives up v andthus the DCO frequency as well. For the design from [20], g2 is quite small, thusv grows gradually until v satisfies fDCO(c,v) = fref . This is the point when phasedifference crosses zero. After that, the trajectory zigzags towards the equilibriumregion. The right sub figure of Figure 3.9 shows the zoom in plot of the small black340.9 0.95 1 1.05 1.12.52.62.72.82.933.13.23.33.4cv0.9 0.905 0.91 0.9152.692.72.712.722.732.74cvFigure 3.9: Trajectory0.9 0.95 1 1.05 1.12.52.62.72.82.933.13.23.33.43.5cvFigure 3.10: Trajectories35box from the left plot. At this time, the BBPFD causes c to oscillate around v/ frefwhile the linear-phase path provides a damping term that keeps these oscillationsmall. Thus, the trajectory has a large number of small zig-zags as v converges toroughly fref ccenter and c tracks v to keep fDCO(c,v)≈ fref and converges to ccenter.Similarly, the red trajectory starts from a place where fDCO is initially greater thanfref , c heads to saturation at cmax; v decreases to fref cmax; and then the trajectoryzigzags to the final locking point. Figure 3.10 shows more trajectories in v× cspace.fDCO(c,v) is modeled by Equation 3.4. fDCO increases with the value of v anddecreases with the value of c. As shown by Equation 3.6, θ∆ changes at a rategiven by the difference between fDCO(c,v) and fref . θ∆ is output by the PFD as thedifference between the time of a rising edge on up and the corresponding risingedge of dn. If the DCO lags the reference, then up will be asserted before dnand the duration of this condition gives the phase difference (negative) between theDCO and reference. Conversely, if the DCO phase is ahead of the reference, dnwill be asserted before up and the duration of this different (now positive) givesthe phase difference between the DCO and reference. The up and dn signals areinput to the DCO as shown in Figure 3.6. When up is asserted without dn, an extrainverter contributes to the drive of each oscillator stage, causing the oscillator torun faster. Conversely, when dn is asserted without up, the oscillator runs slower.This change in the speed of the oscillator translates a duration of up without dn(or vice-versa) on the output of the PFD to a phase shift of the DCO. The ratio ofphase difference output by the PFD to the phase shift of the DCO is the “time gainof the linear phase path and is denoted by K in Equation 3.63.2.3 PFD GlitchHaving developed a simple model for the digital PLL, we considered whether ornot it faithfully captures the behaviour of the real circuit. We realized that thereis a potential “glitch” in the PFD circuit due to the delays involved when resettingthe flip-flops in the PFD circuit. We considered adding these delays to the hybridautomaton model, but doing so would make the model considerably more compli-cated. Our simulations showed that it was hard or even impossible to elicit these36missed Φtrefeventreset θ =  - 1.3 π ∆ θ =  - 1.5π ∆ θ =  - 1.7 π ∆ θ =  + 0.1 π ∆ θ =  - 0.1 π ∆Figure 3.11: PFD Glitchglitches. This prompted a more thorough analysis of the behaviour of the digitalPLL, and I found that I could prove that these PFD glitches are precluded for prac-tical values of the PLL parameters. This allows us to use a simple model of thePFD in the remaining chapters of this thesis where I present my formal verificationof global convergence for the digital PLL. This section examines how delays inthe PFD circuit could lead to glitches and gives constraints on values of the PLLparameters that ensure that such glitches do not occur. An informal explanation forthese constraints is given here, and Appendix A gives the formal proof.Figure 3.11 shows an example of a PFD glitch when the DCO period is slightlylarger than that of the reference. For simplicity, this figure ignores the effect of thelinear phase path (equivalently, it corresponds to a case where the time gain, K, isvery small). The shaded region with each reset pulse shows the time interval duringwhich the flip-flops of the PFD will fail to respond to a rising edge of their clocksdue to the assertion of reset. When the rising edge of ΦDCO occurs very slightlybefore the rising edge of Φref (in other words, θ∆ is close to −2pi), the reset isreleased too late for the PFD to respond to the next rising edge of Φref. Instead, itfirst sees the later rising edge of ΦDCO and then the edge of Φref that occurs slightlylater. For that cycle, dn is asserted before up and the BBPFD will report that theDCO edge occurred before the reference edge. This will cause c to increase, whichincreases the period of the DCO, even though it was already greater than the periodof the reference. We refer to this situation where an edge of one of the clocks is37suppressed by the reset operation as a “PFD-glitch”.To state our condition that ensures there are no PFD glitches, we define a quan-tity, δ which corresponds to the fraction of a period of fref that is not overlappedby PFD reset operations. Let treset be an upper bound for the delay from a risingedge of ΦDCO (or Φref ) that triggers a PFD reset until the reset is complete. Then(1/ fref )− treset is the time during a period of Φref that is not overlapped by resetoperations. Letδ = 1− fref treset (3.7)δ is the fraction of a period of fref that is not overlapped by PFD reset operations.We now state our condition that guarantees that PFD glitches will not occur.Theorem 1. The PFD glitch will not happen if v is in the interval (V blow,V bhi)where:V blow =g2 fref (cmin−ccenter)K + fref cmin− cminKδ ,V bhi =g2 fref (cmax−ccenter)K + fref cmax + cmaxKδ ,The main idea behind these bounds is that, when v is in a fairly wide interval,the linear phase path ensures that θ∆ never crosses the reset threshold δ . Themathematical proof is given in Appendix A. As we can see from the later part ofthe proof, when the linear path parameter K is as big as 0.8 (the value given bythe authors of [20]), the PFD glitch is excluded for most of the parameter valuesfrom Table 3.1. For most of the parameters, V blow is negative which falls out of thelower bound of the valid voltage bound, V bhi exceeds the upper bound of the validvoltage bound. Table 3.2 shows the above bounds for the four set of parameters.Though the set 1’s bounds for v are pretty small, Appendix A shows that it will stillenter into the “safe region”.Table 3.2: Four Sets of v Bounds for PFD Glitchparameters Set 1. Set 2. Set 3. Set 4.v [1.3540, 2.7460] [-0.9, 7.7] [-2.7, 5.5] [-1.3460, 6.0460]38Any verification is only as sound as the model that is used, and that mathe-matical models can only be an approximation of physical reality. For the aboveverification, we modelled the digital PLL as a set of differential equation, while inreality, the c and v integrator are digital components controlled by clocks. There-fore, a time-discretization term is still missing. Because at this point, any changesto the models used in this thesis require lots of manual effort to check all of theresults, we would like to put it as future work currently.3.3 SummaryPhase Locked Loops (PLLs) are common in AMS designs. The key function of aPLL is to adjust the frequency and phase of a voltage-controlled oscillator (VCO)or digitally controlled oscillator (DCO) to track those of a reference signal. Com-mon variations of this design allow PLLs to be used for frequency multiplications,clock-data recovery, phase and frequency modulation and demodulation, and manyother applications. The PLL that we verify in this thesis is a real design fromCICC2010 [20]. It consists of a digitally controlled oscillator (DCO), a phase-frequency detector (PFD), two accumulator, a delta sigma modulator, a low-passfilter (LPF), and a frequency divider between the DCO and the PFD. For simplicity,we omit the divider in our models. This chapter has described each of these compo-nents and presented simple but realistic models for each. We observed conditionsunder which timing glitches could occur in the PFD. We presented conditions thatensure that such glitches never occur, and present a proof in Appendix A that theseconditions are sufficient. The resulting models are the basis for the two approachesto verifying global convergence that we present in the next two chapters.39Chapter 4Verification of GlobalConvergence using SPACEEXSPACEEX is a reachability tool for hybrid system verification. SPACEEX uses sup-port functions to represent continuous states and Linear Hybrid Automata (LHA)models for state computation. This chapter will present the verification of globalconvergence for the digital PLL from Chapter 3 using SPACEEX. Section 4.1presents an ideal SPACEEX model and the over-approximation explosion prob-lem when SPACEEX is used with this ideal model. Section 4.2 presents a revisedapproach that breaks the verification into several separate phases, each with a sim-plified model and works with the dynamics of a real-world problem. The details ofthe two models are included in this section.4.1 Ideal SPACEEX Model and ExplosionSPACEEX uses support functions to represent reachable sets, and dynamics aremodelled by Linear Hybrid Automata (LHA). The linear dynamics can be handledby LHA directly. For nonlinear dynamics, manual work is needed to discretizethe nonlinear part, so that the whole model is Piece-Wise Linear (PWL). The er-ror caused by discretization is then included as differential inclusions part of themodel.The hybrid automata modelling language is called SX in SPACEEX, an XML40like language. Each location is defined by its own differential equations and in-variants. Computation in one location ends on detection of an invariants violation.Transitions are defined by source location, target location, guard condition and as-signment. One location may transit to several locations depending on the guardcondition. The assignment is used when variable reset is needed.4.1.1 Ideal SPACEEX ModelSPACEEX was designed to handle piecewise linear systems and linear inclusions.Most of the digital PLL components have linear models and the non-linearity ofDCO frequency ( fDCO = v/c) appears to be a good candidate for hybridization.Likewise, the steps for v and c are small, which means continuization [3, 4] wouldbe a promising way to handle the discrete time dynamics. SPACEEX had beendemonstrated on a 28-dimensional helicopter control algorithm, which appears tobe scalable. In addition, the user interface is well developed and reasonably welldocumented. Based on these considerations, SPACEEX appears to be a promisingtool for verifying convergence of the digital PLL.Using the SPACEEX modelling language, it would be straightforward to dis-cretize the nonlinear part and model the digital PLL as a PWL hybrid automata.Recall that the behavior model we’ve discussed in Section 3.2.2 is:satc = ((c = cmin)∧ (θδ < 0))∨((c = cmax)∧ (θδ > 0))satv = ((v = vmin)∧ (c > ccenter))∨((v = vmax)∧ (c < ccenter))Resetθ∆ = ((θ∆ ≤−2pi)∧ (θ∆ ≥ 2pi)c˙ = −sign(θ∆)g1 fref , if ¬satc= 0, if satcv˙ = g2 fref (c− ccenter), if ¬satv= 0, if satvθ˙∆ = fDCO(c,v)− fref −Kθ∆ if ¬Resetθ∆reset θ∆ to θ∆−2pisign(θ∆), if Resetθ∆(3.6)The nonlinearity comes from the VCO, that is fDCO. We can discretize fDCO41sat lowNegPossat high!! ≤ !!"#! !! ≥ !!"#!!! > 0! !! < 0!!! < 0!!! > 0!Hybrid automaton for c Integrator! ∆!!! ∆!!! ∆!!! ∆!!!!!!"#$!"# ! !"#$!! !! ≥ !!"#!!! < !!"# ! ! > !!!!!! ≤ !!! !Hybrid automaton for v Integrator!!"# ! !!!! ! !!!! ! !!!! ! !!!! ! !!!! !! < !!!!!! ! < !!!!!! ! < !!!!!! ! < !!!!!!!!!!! !!!!!! > !!!,!! ! > !!!,!! ! > !!!,!! ! > !!!,!! ! > !!!,!! ! > !!!,!!! < !!!,!! ! < !!!,!!!!"# < !!!,! < !!!,! < !!!,! < !!!,! < !!!,! < !!!,! < !!!,! < !!"#!Hybrid automaton for θ∆ IntegratorFigure 4.1: Ideal Hybrid Automata for Each Component of the Digital PLLas follows:fDCO ∈ p1v− p2× c+ p3⊕〈E3 fref 〉 (4.1)Parameters and more details of the derivation are presented in Section 4.2.1.According to this equation we can discretize the fDCO from [20] by separatingthe v space into several overlapping strips and include the related error term asdifferential inclusions. For each interval of v, we computed a linear approximationfor fDCO as a function of c and v, bounded error.Now we can build the PWL model by creating hybrid automata product foreach component. The BBPFD consists two locations controlled by the sign of42θ∆. Because of the digital integrator of c, it also have two saturation locations.Similarly, the PFD has two saturation locations caused by the reset. The parameterv in VCO is discretized and divided into seven overlapping intervals.A simple model of the digital PLL could be obtained by creating the productautomaton from each of the component automata described above. The product ofthe automata shown in Figure 4.1 is the ideal SPACEEX model.4.1.2 ExplosionsAccording to the hybrid automata in Figure 4.1, to verify global convergence, itsuffices to show that this product automaton converges to correct phase and fre-quency lock from all initial states. However, this results in a huge number ofmode transitions. The key issue is that while the output of the PFD could havea mode transition for nearly every cycle, the value of v changes quite slowly. Thus,SPACEEX would need to analyse a large number of mode changes before v settles;in practice, this results in a time out.Furthermore, SPACEEX adds an error term in all directions to the reachablespace with each mode transition. Because v changes very little between modetransitions of the PFD, these error terms overwhelm the convergence of v. All thesewill lead to severe explosion problems in both run-time and the over-approximationof the reachable space when using SPACEEX. Figure 4.2 shows the explosion, forthis plot, we run the model from a tiny box. The value of c changes with the signof θ∆, as the transition frequency increases, c explodes.4.2 Global Convergence VerificationTo avoid the limitations of SPACEEX, we created a model for the digital PLL thatconsisted of two sub-models. The main idea of the verification is to separate thewhole model so that we could consider them according to how close it is to thelocking state. Consider an initial condition where v is much lower than requiredfor lock. As is shown in Figure 3.10, the digital PLL reaches lock in three epochs:• c decreases until it reaches its lower bound.• v increases until fDCO = v/c crosses above fref .43Figure 4.2: SPACEEX Simulation of Limit Cycle Oscillation for Variable c• The trajectories follows a zig-zag path of limit-cycles for c and θ∆ while vclimbs (or drops) to drive c to ccenter.If v is initially much larger than fref ccenter, the trajectory follows a similar se-quence with c saturating at cmax; v dropping until fDCO crosses below fref ; and thena zig-zag to the final equilibrium. To verify global convergence with SPACEEX,we use one model for the first two epochs, and a separate model for the last one.This allows us to overcome the approximation errors for v during the large numberof mode transitions in the final, zig-zag path to equilibrium.In model 1, v is much lower than the equilibrium value. Because c < ccenter,v˙ > 0, from which we conclude that v increases. An equivalent construction whenv is much greater than the equilibrium is employed. For this model, the idealSPACEEX model constructed in Section 4.1.1 is enough to show a bound for v and44θ . Figure 3.10 shows the dynamics of this model.For model 2 (the middle region in Figure 3.10), 1N fDCO oscillates in a smallinterval around fref , and c and v follow a zig-zag path to settle at their equilibriumvalues. Along this path, the interval of θ∆ is small, and it frequently changes itssign, causing a large number of mode transitions that would obscure the progressof v in the SPACEEX analysis. Given the fact that interval for both of | 1Nvc − fref |and θ∆ are small, we can use SPACEEX to show that v is strictly increasing.Together, the results from SPACEEX show that all trajectories that start in thevalid region for model 1 eventually enter the valid region for model 2. All trajec-tories that start in the valid region for model 2, eventually converge to the desiredequilibrium point. Because the union of the valid regions for models 1, 2 coversthe entire space for c, v and θ∆, global convergence of the digital PLL is verified.We will present the details of the two models in the following sections.Table 4.1 shows the sets of parameters that we used for the verification in thisdissertation. Set 2 and Set 3 are the parameters values that we derived from [2, 20]by the author of the design so that it is close to the real design. Set 4 is derivedbased on our behavior model, which is very close to the ones recommended by theauthor. To further explore the digital PLL, we also include a set of parameter withextreme values, that is Set 1.Table 4.1: Four Sets of Parametersparameters Set 1. Set 2. Set 3. Set 4.c [0.9, 1.1] [0.9, 1.1] [0.9, 1.1] [0.9, 1.1]v [1.0, 2.5] [0.2, 3.8] [0.1, 1.5] [1.5, 2.5]θ 5.0 5.0 5.0 5.0fref 2.0 3.0 1.0 2.0ccenter 1.0 1.0 1.0 1.0g1 -0.01 -0.003125 -0.003125 -0.01g2 g1/5 g1/5 g1/5 g1/5K 0.1 0.8 0.8 0.745!"#!! !!"#/!!!!"#/!! !"#!!!!"#!!!"#"$!+1!−1!Figure 4.3: Hybrid Automaton for the PFD4.2.1 Model 1loThis model is for the region where v is much lower than the equilibrium value1.In this region, the sign of phase difference does not frequently changes becausefDCO − fref is still large. Therefore, we can use the simple model from Sec-tion 4.1.1.For the c-accumulator, c changes according to the sign of θ∆ as output by theBBPFD, the c integrator will saturate at the two saturation points, cmin and cmax.Thus, the automaton for c consists of four locations, two locations when c is notsaturated, one with θ∆ positive and c increasing, and the other with θ∆ negativeand c decreasing; plus two saturation locations. The invariants of the negative andpositive location are satisfied when value of c is in interval [cmin,cmax]. There aretwo transitions between these two locations when the phase difference change itssign. When the value of c decreases to cmin and θ∆ is still negative-c won’t decreaseanymore once it reaches cmin, the negative location transitions to the low saturationlocation until the sign of phase difference changes. The positive location is similar.The hybrid automaton for the c accumulator is shown in Figure 4.1.The automaton for the PFD contains one location with a self-loop. This self-transition comes from the reset of the PFD. When the phase difference is higherthan the reset value (or lower than the negative counterpart), 2pi radians are eitheradded or subtracted from θ∆ to bring it close to zero as is shown in Equation 3.6.The transition includes an assignment (reset) to θ∆. Figure 4.3 shows the hybridautomata product of PFD.1This is a model for the case when fDCO lags fref , a similar approach applies for Model 1hi whenfDCO leads fref46The hybrid automaton for θ∆ is shown in Figure 4.1. For model 1, we di-vide v into seven overlapping intervals with the dynamics shown in Equation 4.4.Therefore, there are actually seven sets combined with the above hybrid automa-ton. Combined with the four modes for the c-automaton and the 3 modes for the vautomaton, the product automaton for the PLL has a total of 84 locations.The DCO is more complicated because it has nonlinear dynamics. Non-linearityin SPACEEX is handled by PWL models; thus discretization is needed. Accordingto the behavior model of the DCO in Section 3.2.1, Let c1 = c− c0, v1 = v− v0, c0and v0 is the arithmetic mean of c and v respectively where c∈ [clo,chi], v∈ [vlo,vhi](clo is the lower bound of c, chi is the upper bound of c; Likewise, vlo is the lowerbound of v, vhi is the upper bound of v). Appendix B presents our linearization ofthe fDCO. In particular, from Equation B.8, we get:fDCO ∈ A(1+ v˜i− c˜i) fref ⊕〈E3 fref 〉 (4.2)where,A = 1+E2E1 = c˜v˜E2 = c˜22(1−c˜2)E3 = E1 +(1+ v˜+ c˜+2E1)E2c0 =clo+chi2 ,c ∈ [clo,chi]v0 =vlo+vhi2 ,v ∈ [vlo,vhi]c˜i =c−c0c0c˜ = chi−clochi+clov˜i =v−v0v0v˜ = vhi−vlovhi+vlo(4.3)And,〈x〉 denotes the interval [−x,+x] and x is a scalar; 〈x〉 denotes the enclosing hyper-rectangle for the line segment (−x,+x) if x is a vector.The above equation can be further simplified as follows:fDCO ∈ A(1+ vv0 −cc0) fref ⊕〈E3 fref 〉= A frefv0 v−A frefc0c+A fref ⊕〈E3 fref 〉= p1v− p2c+ p3⊕〈E3 fref 〉(4.4)47Where, p1 =A frefv0, p2 =A frefc0, p3 = A fref .The automata is coded in SPACEEX using SX, and we run the model 1 withthe parameter sets (Table 4.1) in each of the seven regions. We stop the modelwhen c is about to leave the saturation wall, after which it will zig-zag slowly tothe equilibrium. SPACEEX simulation shows a bound for v that is [vlo,vhi] (closeto [ fref cmin, fref cmax]), and a very small bound for θ∆ centred at 0.4.2.2 Model 2loFor model 2, because fDCO is close to fref , θ∆ frequently changes its sign, causinga large number of mode transitions that would obscure the progress of v in theSPACEEX analysis.2 However, we are able to show that the value of fDCO keepsmaking progress. Let w be the difference of the fDCO and fref , that is w = fDCO−fref . Model 2 is constructed as follows:w˙ ∈ fref(g2 +sign(θ∆)g1chmeanw⊕ I1⊕ I2)(4.5)where,I1 = (g1sign(θ∆) fref −g2ccenter)[1/clo,1/chi], if θ∆ > 0= (g1sign(θ∆) fref −g2ccenter)[1/chi,1/clo], if θ∆ ≤ 0I2 = 〈g1ηcwmax〉chmean =2chiclochi + cloηc =chi− clo2chicloThe detailed derivation of above model is presented in Appendix C. The vbound and θ∆ from model 1, the Err term in Equation 4.4 as well as the− fref termin the definition of w are included in the initial value of w. Thus, w0 is presentedas follows:w0 ∈ p1v− p2c+A fref − fref ⊕〈E3 fref 〉 (4.6)2This is a model for the case when fDCO lags fref , a similar approach applies for Model 2hi whenfDCO leads fref48Figure 4.4: Convergence of vNote that for this model, p1 and p2 depend on sign of θ∆.Model 2 has two locations, a negative one and positive one depending on thesign of the θ∆. SPACEEX readily shows that w and θ∆ both converge to intervalsaround 0. Now note that if |w| is small, given v, we can derive a tight bounds forc. This allows us to construct a model using a small interval for w, which showsthat v (and therefore c) converges to its equilibrium value. This model is shown inEquation 4.7, Let Iw be the interval obtained from Equation 4.5. Figure ?? showsthe convergence of v for one set of the parameter.49w = fDCO− fref ∈ p1v− p2c+ p3⊕〈E3 fref 〉− fref = Iwc ∈p1v+ p3− fref − Iw⊕〈E3 fref 〉p2v˙ = g2 fref (c− ccenter)∈ g2 fref(p1v+ p3− fref − Iw⊕〈E3 fref 〉p2− ccenter)(4.7)In Equation 4.7, we compute bounds for c based on the value of w. From these, we,we obtain the inclusion for v˙. As noted earlier, p1, p2, and p3 depend on sign(θ∆);thus, this is a hybrid-automaton model.Unlike the simple model, this modified model does not produce an explosionof the reachable space when analyzed using SPACEEX.Together, these results show that all trajectories that start in the valid regionswill eventually enter the valid region for model 2. Then, all trajectories in model2 converge to the desired small equilibrium region. The regions for model 1 andmodel 2 have covered the entire space for c, v and θ∆. Therefore, we have shownthat the digital PLL will lock to the desired small region starting from anywhere inthe valid space for c, v and θ∆.4.3 SummaryThe ideal SPACEEX model for the digital PLL would be a product automaton ofeach components. My experiments showed that the model caused severe explosionin runtime and over-approximation error due to the large number of transitions.Therefore, we divided the model into two separate models depending on how farthe trajectory is to the desired equilibrium region. Model 1 is in the region wherethe digital PLL is far away from lock. A typical SPACEEX model is built and rununtil the value of c leaves the saturation wall where it is about to zig-zag slowlytowards equilibrium. A bound for v is obtained from this model. Starting fromthis bound, we proved the global convergence by showing that value of fDCO− frefkept decreasing and deriving bounds for c and v that show that they converge to asmall region as well.50Chapter 5Verification of GlobalConvergence usingCOHO-REACHThe difficulties encountered when verifying convergence of the digital PLL us-ing SPACEEX motivated us to consider other verification tools. We choose to useCOHO-REACH [68] because it can compute tighter bounds on the reachable space,and its user interface, a MATLAB API, provides a natural way to decompose theverification problem into a collection of lemmas. This chapter will present ourverification of the global convergence of the digital PLL by using COHO-REACH.Section 5.1 motivates the use of COHO-REACH as well as the its advantages whenused as a lemma prover. Section 5.2 presents the whole proof. More specifically,the verification strategy is first presented. Then in each lemma, the COHO-REACHmodel as well as the verification is described. Section 5.3 concludes this chapter.5.1 Introduction and MotivationAs presented in the last chapter, the digital Phase Locked Loop can be verifiedby SPACEEX given a suitable combination of models. However, the digital PLLis verified in separated phases by different SPACEEX models, which is not wellsupported by SPACEEX. This is mainly caused by SPACEEX’s support function510.96 0.98 1 1.02 1.04 1.06 1.08 1.1 1.12−0.35−0.3−0.25−0.2−0.15−0.1−0.0500.050.10.15cph0.96 0.98 1 1.02 1.04 1.06 1.08 1.1 1.12−0.35−0.3−0.25−0.2−0.15−0.1−0.0500.050.10.15cph(a) (b)Figure 5.1: MATLAB Simulation of One Step of the Limit Cycle Oscilla-tion, (a) Expected, (b) Support Function with Four Directionsrepresentation. The default setting of the support function uses only a few direc-tions, which produce too much over-approximation and the time for the computa-tion grows with the number of directions used.For example, Figure 5.1 shows a single step of a simple model with no modetransition; the plots are from an exact (to within floating point error) computationperformed using MATLAB. Starting from an initial box, as is shown in the figure,the value of c (X axis) should shrink when ph (y axis) approaches 0. The leftplot shows what we expect when ph cross the 0 while the right plot shows whathappened when we use support function with four directions. It is not hard to findthat value of c does not shrink in this case.Accordingly, Figure 5.2 shows the SPACEEX simulation with the exactly samemodel and parameters. The right plot is the one with default number of directions,which took 0.087 seconds. To get the more precise bound for c, we tried 512directions as is presented in the left plot of this figure. The value of c finallyshrinks in this case. However, it took 76.756 seconds for the single step. As theconvergence of the digital PLL takes over thousands of such steps, the computationis extremely expensive.To avoid the numerous zig-zags and excessive over-approximation, we verifiedthe digital PLL by dividing it into several sub-models and validated each of themseparately. The way we used SPACEEX in Chapter 4 was in fact more like a lemmaprover. The digital PLL is verified by the combination of a set of lemmas. In52(a) (b)Figure 5.2: SPACEEX Simulation of One Step of the Limit Cycle Oscillation,(a) uses Support Function with 512 Directions, (b) uses Default Settingsin SPACEEXsome cases, a lemma establishes a post-condition for some flow that is used aspre-condition in other lemmas. SPACEEX provides ‘modules’ (or components),which enables the same structured hybrid automata to have different parametersand combine them to obtain a full model. We could create an automaton for c andan automaton for v. However, the coefficients for the linear model for θ˙∆ depend onboth the mode of the c automaton and the v automaton. SPACEEX doesn’t providea way to “compute” coefficients. Thus, we have to expand the product automatonmanually and write the flow equations for each state of the resulting automaton.This is tedious, so we wrote a JAVA API for generating the XML for SPACEEX.In addition, SPACEEX can tell you what regions are reachable, but it’s hard tofigure out how a trajectory got there. If you could isolate the different transitionsbetween modes, then it would be easier to debug a SPACEEX run.In conclusion, if the system being modelled can be described naturally as aSPACEEX style hybrid automaton, and if SPACEEX can bound the reachable spaceby its fix point computation, then, it would be a great tool. For the digital PLL,these conditions weren’t satisfied. We can decompose the verification problem intosub-problems that SPACEEX can handle, but SPACEEX provides no mechanism forcomposing these “lemmas”.53The issues described above motivate us to consider other tools which are moresuitable to use as lemma prover, with tighter error bounds and more convenientto debug. COHO-REACH is an outstanding option, it’s excellent in the two mainfeatures that we are looking for: tight error bounds and convenience to use as alemma prover.COHO-REACH is a reachability tool for verifying dynamic systems. It cananalyze non-linear hybrid automata. The user describes continuous flows by pro-viding a “model function”. This function takes as an argument bounds on the partof the phase-space under consideration and returns a linear differential inclusionthat contains the actual dynamics in that region. For non-linear dynamic systems,non-linear terms are linearized and hybridization applies.As a hybrid automata abstraction, COHO-REACH provides modes and switch-ing surfaces. Modes are defined by its dynamics and invariants. Computation forone mode ends on detection of a violation of an invariant. Transitions are thenchecked to find the next computation mode.COHO-REACH allows the user to define termination conditions for the analysis.While COHO-REACH can compute fix points in a manner similar to SPACEEX,users are able to use the feature of user-defined termination conditions. This allowsone to verify a collection of “pre-condition→ post-condition” style lemmas, thattaken together prove the main result. This is helpful because users are able tostate these pre- and post-conditions directly in Matlab code. COHO-REACH alsoprovides functions that make it easy for the user to check that the reachable regionat each step, or at a termination condition or transition, satisfies a predicate suchas containment in a hyperrectagle or projectagon (see Section 2.2.2). Users don’tneed to manually inspect reachable regions to determine whether or not the modelhas the desired properties as a user of SPACEEX must.Unlike SPACEEX, COHO-REACH performs on-the-fly linearization to constructinclusions as they are needed. This not only produces tighter inclusions; it alsoavoids switching events that are merely artifacts of the hybridization process andnot actual switching behaviours of the underlying system. In addition, COHO-REACH determines support vectors by adaptively computing the time-advancednormals of the polyhedron’s faces. Moreover, COHO-REACH uses an error-boundwith direction, and adapts the number of terms in its approximations to achieve54a desired tolerance. In contrast, SPACEEX uses a simple norm-based method tohandle inclusions and finite approximation of matrix exponentials.COHO-REACH is run in MATLAB environment, most of the COHO-REACHcode is written in MATLAB and handy to access. Therefore, users are able tocall any MATLAB built-in functions and use the debugging facilities that MAT-LAB provides. More importantly, on-the-fly linearization frees us from writingrepeated blocks of cut-and-paste code as needed by SpaceEx. Of course, we coulduse our on-the-fly linearization model to generate the inclusions used in SpaceEx,but we would still have to generate the specific inclusions required for the statichybridization, and write XML code for each of these hybridization. Coho-Reachsimply needs the model function and calls it for each linearization that it needs.This provides a very flexible way to produce inclusions from non-linear models.The reasons described above motivate us to use COHO-REACH as a lemmaprover when verifying the global convergence of the digital PLL. In the followingsections, we will show how the COHO-REACH models are built and how COHO-REACH can be used to prove the lemmas.5.2 Global Convergence VerificationAs with SPACEEX, a COHO-REACH model is defined by its modes, invariants andcorresponding transitions. Furthermore, COHO-REACH allows for user-specifiedstop conditions and time step. For each mode, flows are specified by a user-definedmodel as described in Section 5.1. These models take a region as an argument andreturn a linear-differential inclusion that contains all flows in that region. A transi-tion connects two modes, the source mode and the target mode. The transition maybe triggered by the violation of the one or more invariants. A “resetMap” functionis provided to change the initial condition of the target state. COHO-REACH alsoprovides user-defined stop conditions that allow the user to inspect the reachableregion when the condition is reached. We use these features to decompose the ver-ification into a collection of lemmas that are checked by COHO-REACH. Further-more, the logical dependencies of these lemmas are checked by our Matlab scriptsusing functions provided by the COHO-REACH API. As noted above, there are log-ical dependencies between these lemmas; for example, one lemma may show that55all trajectories eventually enter some region, and another lemma has a hypothesisthat all trajectories have entered that region. We explicitly test these relationshipsin our Matlab scripts using functions provided by the Coho-Reach API.We verified global convergence of the digital PLL using COHO-REACH byproving the following four lemmas:Lemma 1: Starting anywhere in the valid initial space (i.e. c, v, and θ∆ boundedby their hard limits), trajectories either reach a point where (c= cmin)∧(θ∆ <0) or (c = cmax)∧ (θ∆ > 0) or θ∆ = 0.Lemma 2: In saturation modes (modes with (c = cmin)∧ (θ∆ < 0) and (c =cmax)∧ (θ∆ > 0)), all the trajectories will finally leave the mode with θ∆small.Lemma 3: All trajectories that cross the θ∆ = 0 plane eventually enter a narrowstrip parallel to the plane vc = fref .Lemma 4: For trajectories that are close enough to the plane vc = fref , each zig orzag brings it closer to a small region that contains the final limit cycle withc≈ ccenter.Together, these four lemmas show that starting from any initial point, the digitalPLL eventually reaches proper lock.5.2.1 Initial TrajectoriesLemma one: Starting anywhere from the valid space, trajectories either reach apoint where (c = cmin)∧ (θ∆ < 0) or (c = cmax)∧ (θ∆ > 0) or θ∆ = 0.Recall that the behavior model we’ve discussed in Section 3.2.2 is:56PosNegSat High Sat Lowtermination  conditiontermination  conditionFigure 5.3: COHO-REACH Modes and Switching Surfacessatc = ((c = cmin)∧ (θδ < 0))∨((c = cmax)∧ (θδ > 0))satv = ((v = vmin)∧ (c > ccenter))∨((v = vmax)∧ (c < ccenter))Resetθ∆ = ((θ∆ ≤−2pi)∧ (θ∆ ≥ 2pi)c˙ = −sign(θ∆)g1 fref , if ¬satc= 0, if satcv˙ = g2 fref (c− ccenter), if ¬satv= 0, if satvθ˙∆ = fDCO(c,v)− fref −Kθ∆ if ¬Resetθ∆reset θ∆ to θ∆−2pisign(θ∆), if Resetθ∆(3.6)These equations can be linearized (see Equation 4.4). The reachable regionis recorded as it crosses a switching surface to further modes and used to derivethe pre-condition part of subsequent lemmas. As is shown in Figure 5.3, modescomputations are ended on user-defined termination conditions, those switchingsurfaces are connected to other reachable modes.Algorithm 1 and Algorithm 2 shows the pseudo-code of the main functionsused for verifying this lemma. Function Run mode1() defines the initial projec-57tagon initPh, termination condition callBacks.exitCond, and switching surfacesph sat low, ph sat high, ph mode2. Line 1−5 initiate the projectagon by bound-ing box bbox. Line 9− 10 set the termination conditions and invariants. Line 11calls the Function mode1() from Algorithm 2 and returns the reachability data andcomputation time. If computation time exceeds the maximum computation timemaxCompT , then, it will stop and return False as Line 13− 16 specified. Other-wise, we go through 16−26, get the post condition PostCond from switching sur-faces. Data are collected by COHO-REACH function ph union() and ph intersect().Function mode1() calls COHO-REACH function ha reach() to compute thereachable projectagon given the initial condition, invariants, termination conditionand parameters. The hybrid automaton (States and transitions) is defined in Func-tion EX PLL MODE1 and the real model (Differential equations) are described inFunction PLL MODE1 from Line 12−16 in Algorithm 2.COHO-REACH readily validates this lemma. See Appendix D for more detailsof the code.5.2.2 Leaving the WallLemma two: In saturation modes (modes with (c = cmin)∧ (θ∆ < 0) and (c =cmax)∧ (θ∆ > 0)), all the trajectories will finally leave the mode with θ∆ small.The COHO-REACH model includes the two saturation modes, where c˙ = 0,because the value of c is saturated at the cmin or cmax, the nonlinear term disap-pears and the model is purely linear. In this case, we can treat c as a parameter(fixed) rather than a variable (changing) in the model. The hybrid automaton canbe constructed according to the following ordinary differential equations:v˙ = g2 fref (csat − ccenter)θ˙∆ =vcsat− fref −Kθ∆(5.1)csat can be cmin or cmax depending on the saturation mode. For this model, wedetermine an interval for v as trajectories leave the wall, i.e. when θ˙∆ reaches 0.When a trajectory leaves the wall, it zig-zags towards the equilibrium state. Theinterval for v along with the values for c and θ∆ should satisfy the precondition for58Algorithm 1 Run mode1Input: initPh, pars // initPh:initial projectagon, pars: parametersOutput: Pass (boolean), PostCond1: if nargin < 1 then2: bbox = [[0.9,1.1]; [0.2,3.8]; [−2pi,+2pi]] //bounds for c,v and θ∆3: initPh = ph createByBox(dim, planes,bbox) //create initial...projectahedron4: initPh= ph convert(initPh, ‘convex′) //mark the reachable region as convex,...COHO-REACH will optimize for this...case5: end if6: if nargin < 2 then7: Pars = getPLLParameters() //use default PLL parameters8: end if9: callBacks.exitCond = @(info)ex mode1 exitCond(in f o, pars.maxCompT )//terminate reachability computation on time-out10: invs=@ex mode1 invariants() //set mode invariants: cmin < c < cmax11: ha=mode1(callbacks, initPh, invs, pars) // create the hybrid automaton note, it...has one mode, transitions to the...other modes are handled by user...code12: [data, time] = LoadData(ha)13: if time≥ pars.maxCompT then14: Pass←False //time out15: PostCond← [ ]16: else17: Pass←True // computation completed, compute initial regions for other...checks18: ph sat low = ph createByBox(dim, planes, [...]) // initial region for...c = cmin19: ph sat high = ph createByBox(dim, planes, [...]) // initial region for...c = cmax20: ph mode2= ph createByBox(dim, planes, [...]) // initial region for θ∆...crossing 021: for i= 1 : length(data) do // gather regions for entering other modes from...reachable region for each time step22: sat low = ph union(sat low, ph intersect(data(i), ph sat low))23: sat high = ph union(sat high, ph intersect(data(i), ph sat high))24: mode2 = ph union(mode2, ph intersect(data(i),mode2))25: end for26: PostCond← [sat low; sat high; mode2]27: end if59Algorithm 2 Mode1Input: callBacks, initPh, pars //initPh came from mode 1Output: ha //reachData1: ha = EX PLL MODE1(initPh,callBacks, invs, pars)2: ha = ha reach(ha)3: function EX PLL MODE1(initPh,callBacks, invs, pars)4: phOpt = @getSettings()5: States(1)= ha state(‘name′,@(l p)(PLL MODE1(pars, l p)), invs, phOpt,6: callBacks)7: States(2) = ...8: trans = @getTransitions()9: Source = ‘name′10: ha = ha create(‘name′,states, trans,scource, initPh)11: return ha12: end function13: function PLL MODE1(pars, lp)14: A = [...] // define the differential inclusions for the dynamics systems,...linearize it if needed. See Appendix B for more details of the...linearization15: B = [...]16: U = [...]17: return int create(A,B,U)18: end functionLemma 3.This lemma can be verified by COHO-REACH in a similar fashion as the pre-vious lemma. This model only needs two variables, the number of dimensions andprojection planes are different than used to prove lemma 1. Again, COHO-REACHreadily proves this lemma. See Appendix D for the detailed code.5.2.3 Convergence of fDCO to frefLemma three: All trajectories that cross the θ∆ = 0 plane eventually enter a narrowstrip that parallel to the plane vc = fref .Trajectories in which c saturates were addressed by Lemma 2; therefore, wecan exclude saturation here and model the BBPFD using only two modes (modesNeg and mode Pos from Figure 5.3). The trajectories under consideration repeat-60edly cross the θ∆ = 0 plane in alternating directions. At each such crossing, weconsider |v− fref c| as a measure of the “distance” of the trajectory from the idealvalue. We show that this distance decreases rapidly with successive crossings untilthe trajectory is in a narrow strip entered along v/c = fref .Let ds(c,v) be the distance of (c,v) to the line v/c = fref ; more precisely,ds(c,v) = v− c fref . To initialize the model, rectangles that are parallel to the vc =fref plane are constructed. COHO-REACH shows that ds will keep shrinking untilds ∈ 〈0.03〉 from the model parameters from set 1, Table 3.1.Algorithm 3 shows the pseudo-code for discharging this lemma. We knowthat, for trajectories below vc = fref , v ≥ (c−ds) fref ; thus, the maximum distancedslo can be calculated as shown on Line 7 of Algorithm 3. On Line 8, stripNumdefines the number of strips. Each strip corresponds to an interval for ds. Line 9defines an array stripBoundLo where stripBoundLo(i) gives the lower bound fords for strip i. stripBoundLo equally slice it into stripNum strips. For each suchstrip, Line 11−12 construct a projectagon from the constraints of a linear programthat describe the strip. Unlike the previous termination condition, we have COHO-REACH compute the reachable set until the projectagon for this mode is empty(i.e. all trajectories have transitions to the final mode described in Section 5.2.4).Given the data, Line 21 checks whether the distance decreases. The lemma isdischarged if and only if ds decreases by at least the prescribed amount for eachstrip. By doing this we can show that the projectagon keeps moving upwards alongthe plane vc = fref . Likewise, we can construct the higher part. See Appendix D formore details of the code.5.2.4 Convergence of c to ccenterLemma four: For trajectories that are close enough to the plane vc = fref , each zigor zag brings it closer to a small region that contains the final limit cycle withc≈ ccenter.The basic idea of this model is similar to our second model of SPACEEX ver-ification in Section 4.2.2 but with model specifically designed for COHO-REACH.Because θ∆ frequently changes its sign, causing a large number of mode transitionsthat would obscure the progress of v, we show that starting from the narrow strip,61Algorithm 3 Lemma3Input: ds, parsOutput: Pass(boolean)1: if nargin < 1 then2: ds = 0.03 // defined default value for ds3: end if4: if nargin < 2 then5: Pars = getParameters()6: end if7: dslo = (cmax fref − vmin)/ fref8: stripNum = 5;9: stripBoundLo=(ds(stripNum :−1 : 0)+dslo(0 :+1 : stripNum))./stripNum10: for k = stripNum :−1 : 1 do11: ds1 = stripBoundLo(k+1)12: ds2 = stripBoundLo(k)13: init low = l p create([1,−1/ fref ,0],ds1)14: init hi = l p create([−1,1/ fref ,0],−ds2)15: init = l p and(init low, init hi)16: initPh= ph createByLP(dim, planes, l p and(init, l p createByBox(bbox)))17: callBacks.exitCond =@(in f o)ex mode1 exitCond(in f o, pars.maxCompT )18: invs = @ex mode1 invariants()19: ha = mode1(callbacks, initPh, invs, pars)20: data = LoadData(ha)21: Pass←@(data)isDsDecreasing(data)22: end foreach mode transition will bring it closer to the equilibrium point compared with thelast step. The model can therefore be constructed by shifting both value of v and cso that they are centered at vc = fref . By doing this, it becomes straightforward tosee the change of c and v for each transition.v1 = v− vmid c1 = c− cmidvmid = 12(vlo + vhi) cmid =12(clo + chi)(5.2)The model can therefore be formulated as follows:62c˙1 =−sign(θ∆)g1 frefv˙1 ∈ g2 fref c1−g2 fref (ccenter−vmidfref)⊕g2vmid v˜θ˙∆ ∈ A(1+v1vmid−c1cmid) fref − fref −Kθ∆⊕E3 frefA = 1+E2E1 = c˜v˜E2 =c˜22(1− c˜2)E3 = E1 +(1+ v˜+ c˜+2E1)E2c˜ =chi− clochi + clov˜ =vhi− vlovhi + vlo(5.3)[clo,chi] and [vlo,vhi] are the intervals for c and v respectively. The derivation ofthis model is formally proved by Z3, which is presented in Appendix B. Ourexperience was that it is common to make minor, algebraic errors in the derivationsof the inclusions. By verifying the results with an SMT solver, we ensure thatthe inclusions over-approximate the behaviors of the real system as intended. Inprinciple, it would be good to verify all of the models presented in this thesis inthis manner. Our use of Z3 for this purpose was introduced late in the research. Wepresent it as a proof-of-concept result, and note that such verification of manuallyderived formulas has much wider application that just checking this one formula.For the model described above, our ideas is to consider a narrow parallelogramin c× v space whose long edges are parallel to the vc = fref line. From simulation,we expect c to zig-zag in this region while v moves slowly towards fref ccenter andc tracks v and converges to ccenter. Each “zig” or “zag” turns around when θ∆crosses 0. We identify a bound on the magnitude of θ∆ and | fDCO(v,c)− fref | atthese turn-around points. We show that at each successive turn-around point thesebounds are satisfied and that v moves closer to the target value of fref ccenter. The63key observation is that we can write formulate a linear inclusion that applies to theentire parallelogram or large pieces of it, and then re-write the dynamics to applyto a single zig or zag of the trajectory. By showing that every zig or zag of thetrajectory makes progress towards the equilibrium region, we prove convergence.In more detail, let vˆ be the value of v when θ∆ changes sign. If vlo and vhiare bounds for v in the parallelogram, then vlo ≤ vˆ ≤ vhi. Let cˆ = vˆ/ fref . Let c(t)and v(t) be the values of c and v at t time units after turn-around point. We define∆c(t) = c(t)− cˆ and ∆v(t) = v(t)− vˆ. We can rewrite the dynamics of the PLLin terms of a linear differential inclusion for ∆c, ∆v and θ∆, where the inclusionis valid for any vlo ≤ vˆ ≤ vhi. We now use COHO-REACH to analyse a trajectorythat starts with a ∆c(0) and θ∆(0) in the bounds identified above. COHO-REACHcomputes the reachable region at the next sign change of ∆θ plus some uncertaintytime to model that the phase detector may not detect sign change for up to onecycle of fref . We include a check that the values for ∆c and θ∆ are within theprescribed bounds at this turn-around point and that ∆v has made progress towardsthe equilibrium value.For this model COHO-REACH runs several intervals for v in the narrow stripfrom the previous step for both mode 1 and mode 2. The result shows that allthe trajectories in the narrow strip will finally converges into a small box, whichcontains the exact equilibrium point.Together, these lemmas show that all trajectories will go to the wall when csaturates or move towards the plane with vc = fref and a small interval for θ . Thosetrajectories for which c saturates will eventually leave the wall and enter into anarrow strip centered about v/c = fref and a small interval for θ . Finally, startingfrom the narrow strip and the small interval for θ , COHO-REACH is able to showthat each zig or zag will bring it closer to the equilibrium state. Combining all theselemmas, we therefore, claim that the digital PLL will lock to the desired frequencyand phase with c close to ccenter starting anywhere from the valid space. Currently,we check this composition of lemmas manually. Automating this process is a topicfor future research.645.3 SummaryWe choose COHO-REACH as for verifying the digital PLL because it providestighter error bounds than SPACEEX, and its programmable interface allowed us tonaturally decompose the verification into lemmas that could be verified by reacha-bility computation. In this chapter, we divide the verification work into four lem-mas. Together, these lemmas show that all trajectories will go to the wall when csaturates or move towards the plane with vc = fref and a small interval for θ . Thosetrajectories for which c saturates will eventually leave the wall and enter into anarrow strip centered about v/c = fref and a small interval for θ . Finally, startingfrom the narrow strip and the small interval for θ , COHO-REACH is able to showthat each zig or zag will bring it closer to the equilibrium state. Combining theselemmas, we have shown that starting from anywhere in the valid initial space, thedigital PLL will lock to the desired frequency and phase with c close to ccenter.65Chapter 6Conclusion and Future WorkThis dissertation presents the global convergence verification of a digital phaselocked loop using two leading reachability computation tools: SPACEEX and COHO-REACH. Each component of the digital PLL is modelled by a hybrid automata,and the global convergence property is verified using reachability techniques. Sec-tion 6.1.1 summarizes the SPACEEX verification; and Section 6.1.2 concludes theCOHO-REACH verification. Section 6.1.3 gives a brief of these two tools usedin terms of error bounds, usability, implementation etc. Section 6.2 presents thefuture work.6.1 ConclusionA phase locked loop is a typical AMS circuit; it is a ubiquitous component in SoCdesigns and widely used in analog and mixed signal designs. The design of PLLsare of great importance because a failure to lock would block further test of entirechip or major subsystem. Today, device scaling to smaller feature sizes motivatesa shift from traditional, analog designs to “digital” architectures. The PLL verifiedin this dissertation is a state-of-the-art digital one from [20].The digital PLL verified in this thesis consists of a digitally controlled oscillator(DCO), a PFD, a bang-bang phase frequency detector (BBPFD), low pass filter(LPF) and a frequency divider between DCO and PFD. We modelled the PLL asa set of differential equations. To better capture the behavior of the real circuit,66we’ve carefully analyzed a “glitch” in the PFD circuit caused by the delay duringresetting of the flip-flops. It turns out that for a wide interval of the control voltage,v, the linear phase path ensures that this PFD glitch never occurs. We presentedconditions and a proof in Appendix A. The behavior model that we derived is usedas the basis for verifying global convergence.Hybrid automata provide an appealing abstraction for the modelling and verifi-cation of control systems and analog-mixed-signal circuits. SPACEEX and COHO-REACH are two leading verification tools implementing this abstraction but theydiffer in many ways. This dissertation presented the global convergence verifica-tion of the digital PLL first using SPACEEX and then using COHO-REACH.6.1.1 Global Converge Verification by SPACEEXSPACEEX is designed from the outset to handle large linear hybrid automaton mod-els. The piecewise linear inclusions model the PLL components quite well. Ac-cording to the behavior model, the main nonlinear term comes from the DCO, forwhich hybridization applies. Therefore, an ideal SPACEEX model for the digitalPLL would be an automata product of each components.However, we encountered problems with long compute times and large over ap-proximations when SPACEEX computed non-trivial fix points for cycles of modes.Instead, we organize the modes of the automaton according to the typical behaviorof the PLL during lock to avoid cycles of modes. With this change, the globalconvergence of the digital PLL is verified.6.1.2 Global Converge Verification by COHO-REACHMotivated by the need for a reachability tool with tighter bounds and flexibility, wethen verified the PLL using COHO-REACH. COHO-REACH was designed with afocus on algorithms for efficiently computing tight bounds on the reachable regionsfor weakly non-linear hybrid automata.We showed how differential inclusions for modelling the PLL can be obtainedby linearizing the behavior model. Because COHO-REACH uses on-the-fly lin-earization, the error term of the inclusion are specific for the current time-step,and smaller than those for a fixed-partitioning such as used by SPACEEX. Further-67more, on-the-fly linearization avoids mode transitions that are only an artifact ofhybridization. Global convergence is established for the PLL by proving a set offour lemmas using COHO-REACH. In particular, these lemmas formalize simpleobservations which indicate where the trajectories go from an initial projectagon.From our experience, COHO-REACH shows great potential to be used as a lemmaprover.6.1.3 Pros and Cons of SPACEEX and COHO-REACHSPACEEX and COHO-REACH are both well developed reachability tools for veri-fying control systems. They both implement the hybrid automata abstraction anduse Piece Wise Linear (PWL) model for nonlinear dynamics. However, they havesignificant differences when applied to a realistic verification problem.RepresentationSPACEEX uses support function to represent polyhedra while COHO-REACH useprojectagon to represent polyhedra.For SPACEEX, the default choice for support functions uses only a few di-rections which produces excessive over-approximation. On the other hand, using asufficient number of support function directions to verify the PLL results in imprac-tically large run-times for SPACEEX. Figure 6.1 shows the SPACEEX simulationof a single step of the zig-zag with 512 directions and 2 directions respectively.The right plot only took 0.087 seconds while the left plot took 76.756 seconds.Because verifying the PLL requires thousands of such steps, it would take unac-ceptably long to show that the digital PLL converges.LinearizationBoth SPACEEX and COHO-REACH use PWL models for nonlinear dynamics. How-ever their implementations are different: SPACEEX uses hybridization with prede-termined partitioning of the state-space into linear, differential inclusion models.For each region, non-linearities are captured by an inclusion term that is largeenough to account for the worst-case discrepancy with the linear approximation foreach region. On the other hand, COHO-REACH uses on-the-fly linearization and68(a) (b)Figure 6.1: SPACEEX Simulation, (a) uses Default Settings in SPACEEX, (b)uses Support Function with 512 Directionsadaptive construction of support vectors. On-the-fly linearization, as an alternativeto hybridization for producing linear inclusions from non-linear models, producestighter inclusions by tailoring them to the regions in which they are needed. In ad-dition, by avoiding extra mode switching between different regions, the on-the-flyapproach further reduces over-approximation and produces a more efficient algo-rithm.TransitionsFor modes (locations, states) transitions, SPACEEX and COHO-REACH both followthe basic hybrid automata model, that is, the computation stops on detection ofinvariant violation. COHO-REACH provides more flexibility as it supports user-defined stop conditions. This feature was of great use in our verification. However,SPACEEX handles the transitions more accurately by performing both forward andbackward computation to bound the next starting polyhedra. COHO-REACH onlyperforms the forward computation.UsabilitySPACEEX has a well developed graphic interface and XML-like language for mod-elling while COHO-REACH is run in MATLAB environment. SPACEEX provides69‘modules’ (or components), which enable same structured hybrid automata to havedifferent parameters and integrate them into high order differential equations. How-ever, SPACEEX does not allow the parameters of one component model to be func-tions of the state of other models. This prevented us from using SPACEEX’s com-ponent facility when modelling the digital PLL. Instead, we implemented a simpleJava API to generate flattened, hybrid-automata descriptions in SPACEEX’s XMLlike language from structured Java code. Another limitation of SPACEEX is thatit cannot separate a verification task into separate phases. In particular, it cannotuse the post-condition from one analysis as a pre-condition for further reachabilitycomputation. Instead, such relationship had to be checked manually by inspectingplots generated by SpaceEx.In addition, SPACEEX is hard to debug, it is difficult to follow the flow espe-cially when there are several different transitions from one location.COHO-REACH is run in MATLAB environment, most of the COHO-REACHcode is written in MATLAB and handy to access. Therefore, users are able to callany Matlab built-in functions and the debug functions that Matlab provides. Moreimportantly, on-the-fly linearization frees us from writing chunks of repeated codeand provides us with significant flexibility. Thus, we believe user would benefitmore by having more control over the code.6.2 Future WorkThere are many areas for future work.• We would like to add quantization error to our models. This error arisesfrom approximating discrete accumulators with continuous integrators. Webelieve that such error terms could be automatically derived and added to theCOHO-REACH models. The c and v integrator are both digital parts whichupdates value by the controlled clock, the continuous differential equationsignores this digital feature.• We would like to complete models for low-pass filter and Delta-Sigma mod-ulator. The low-pass filter seemed like an obvious candidate to include in ourmodel: we can model it as a purely linear system with three state variables70with no mode transitions. It would be interesting, because the filter’s cut-offfrequency for the digital PLL is much higher than the corner frequency ofthe accumulators. This results in a stiff model. On the other hand, we couldalso consider the case where the cut-off frequency of the low-pass filter isclose to that of the v-integrator. In this case, the PLL will be unstable, andour verification should show this.• As our COHO-REACH verification uses COHO-REACHas a lemma prover,we would like to develop a lemma prover tool based on COHO-REACH.Chapter 5 shows that most of the lemmas indicate where the trajectories gofrom initial states. Thus, it is quite promising to abstract this pattern.• Another interesting property of PLL verification is the lock time. We wouldlike to provide bounds on lock time (excluding metastability) in the future.71Bibliography[1] G. Al-Sammane, M. Zaki, and S. Tahar. A symbolic methodology for theverification of analog and mixed signal designs. In Design, Automation Testin Europe Conference Exhibition, 2007. DATE ’07, pages 1–6, April 2007.doi:10.1109/DATE.2007.364599. → pages 13[2] E. Alon. Parameter values for the digital PLL. personal communication,2013. → pages iv, 45[3] M. Althoff, S. Yaldiz, A. Rajhans, X. Li, B. Krogh, and L. Pileggi. Formalverification of phase-locked loops using reachability analysis andcontinuization. In Computer-Aided Design (ICCAD), 2011 IEEE/ACMInternational Conference on, pages 659–666, Nov 2011.doi:10.1109/ICCAD.2011.6105400. → pages 6, 20, 41[4] M. Althoff, A. Rajhans, B. H. Krogh, S. Yaldiz, X. Li, and L. T. Pileggi.Formal verification of phase-locked loops using reachability analysis andcontinuization. Commun. ACM, 56(10):97–104, 2013. → pages 6, 20, 41[5] R. Alur. Predicate abstraction for reachability analysis of hybrid systems.ACM Trans. Embedded Comput. Syst, 5:2006, 2006. → pages 15[6] R. Alur and D. L. Dill. A theory of timed automata. Theoretical ComputerScience, 126:183–235, 1994. → pages 14[7] R. Alur, T. Dang, and F. Ivancic. Counter-example guided predicateabstraction of hybrid systems, 2003. → pages 15[8] E. Asarin, T. Dang, and A. Girard. Reachability analysis of nonlinearsystems using conservative approximation. In O. Maler and A. Pnueli,editors, Hybrid Systems: Computation and Control, volume 2623 of LectureNotes in Computer Science, pages 20–35. Springer Berlin Heidelberg, 2003.ISBN 978-3-540-00913-9. doi:10.1007/3-540-36580-X 5. → pages 4, 12,1572[9] E. Asarin, T. Dang, and A. Girard. Reachability analysis of nonlinearsystems using conservative approximation. In In Oded Maler and AmirPnueli, editors, Hybrid Systems: Computation and Control, LNCS 2623,pages 20–35. Springer-Verlag, 2003. → pages 15[10] E. Asarin, T. Dang, and A. Girard. Hybridization methods for the analysis ofnonlinear systems. ACTA INFORMATICA, 43:451–476, 2007. → pages 15,18[11] J. D. B. Dipert and W. Strauss. The design and verification challenge for thenext decade. http://low-powerdesign.com/fosler designandverify.htm. →pages 3[12] R. Bagnara, E. Ricci, E. Zaffanella, and P. Hill. Possibly not closed convexpolyhedra and the Parma polyhedra library. In Static Analysis, volume 2477of Lecture Notes in Computer Science, pages 213–229. Springer BerlinHeidelberg, 2002. ISBN 978-3-540-44235-6.doi:10.1007/3-540-45789-5 17. → pages 16[13] A. Balivada, Y. Hoskote, and J. Abraham. Verification of transient responseof linear analog circuits. In VLSI Test Symposium, 1995. Proceedings., 13thIEEE, pages 42–47, Apr 1995. doi:10.1109/VTEST.1995.512615. → pages10[14] B. Boser and B. Wooley. The design of sigma-delta modulationanalog-to-digital converters. Solid-State Circuits, IEEE Journal of, 23(6):1298–1308, Dec 1988. ISSN 0018-9200. doi:10.1109/4.90025. → pages 2[15] R. Bryant. Graph-based algorithms for boolean function manipulation.Computers, IEEE Transactions on, C-35(8):677–691, Aug 1986. ISSN0018-9340. doi:10.1109/TC.1986.1676819. → pages 10[16] C. Chien. Digital Radio Systems on a Chip: A Systems Approach. KluwerAcademic Publishers, Norwell, MA, USA, 2001. ISBN 0-7923-7260-3. →pages 2[17] A. Chutinan. Hybrid system verification using discrete modelapproximations. PhD thesis, Carnegie Mellon University, 1999. → pages 4,12, 16[18] E. Clarke and E. Emerson. Design and synthesis of synchronizationskeletons using branching time temporal logic. In D. Kozen, editor, Logicsof Programs, volume 131 of Lecture Notes in Computer Science, pages7352–71. Springer Berlin Heidelberg, 1982. ISBN 978-3-540-11212-9.doi:10.1007/BFb0025774. → pages 11[19] P. Cousot and N. Halbwachs. Automatic discovery of linear restraints amongvariables of a program. In Proceedings of the 5th ACM SIGACT-SIGPLANSymposium on Principles of Programming Languages, POPL ’78, pages84–96. ACM, 1978. doi:10.1145/512760.512770. → pages 4[20] J. Crossley, E. Naviasky, and E. Alon. An energy-efficient ring-oscillatordigital PLL. In Custom Integrated Circuits Conference (CICC), 2010 IEEE,pages 1–4, Sept 2010. doi:10.1109/CICC.2010.5617417. → pages iv, 26,27, 33, 34, 38, 39, 42, 45, 66[21] N. Culp. Digital phase locked loop, Mar 1986. US Patent 4,577,163http://www.google.com/patents/US4577163. → pages 2[22] T. Dang. Verification and synthesis of hybrid systems. PhD thesis, InstitutNational Polytechnique de Grenoble, 2000. → pages 4, 12, 15[23] T. Dang, A. Donz, and O. Maler. Verification of analog and mixed-signalcircuits using hybrid system techniques. In A. Hu and A. Martin, editors,Formal Methods in Computer-Aided Design, volume 3312 of Lecture Notesin Computer Science, pages 21–36. Springer Berlin Heidelberg, 2004. ISBN978-3-540-23738-9. doi:10.1007/978-3-540-30494-4 3. → pages 12[24] T. Dang, C. Guernic, and O. Maler. Computing reachable states fornonlinear biological models. In P. Degano and R. Gorrieri, editors,Computational Methods in Systems Biology, volume 5688 of Lecture Notesin Computer Science, pages 126–141. 2009. ISBN 978-3-642-03844-0.doi:10.1007/978-3-642-03845-7 9. → pages 18[25] L. De Moura and N. Bjørner. Z3: An efficient smt solver. In Proceedings ofthe Theory and Practice of Software, 14th International Conference on Toolsand Algorithms for the Construction and Analysis of Systems,TACAS’08/ETAPS’08, pages 337–340, Berlin, Heidelberg, 2008.Springer-Verlag. ISBN 3-540-78799-2, 978-3-540-78799-0. → pages 87[26] I.-S. Dhingra. Formalising an integrated circuit design style in higher orderlogic,. Ph.D. dissertation, Computer Laboratory, Cambridge University,Nov,1988. → pages 19[27] Z. J. Dong, M. Zaki, G. Al-Sammane, S. Tahar, and G. Bois. Checkingproperties of PLL designs using run-time verification. In Microelectronics,742007. ICM 2007. International Conference on, pages 125–128, Dec 2007.doi:10.1109/ICM.2007.4497676. → pages 6, 19[28] A. Eggers, M. Fra¨nzle, and C. Herde. SAT modulo ODE: A direct SATapproach to hybrid systems. In Automated Technology for Verification andAnalysis, volume 5311 of Lecture Notes in Computer Science, pages171–185. Springer Berlin Heidelberg, 2008. ISBN 978-3-540-88386-9.doi:10.1007/978-3-540-88387-6 14. → pages 16[29] M. Fra¨nzle. Analysis of hybrid systems: An ounce of realism can save aninfinity of states. In CSL, volume 1683 of LNCS, pages 126–140. Springer,1999. → pages 16[30] G. Frehse. Compositional verification of hybrid systems using simulationrelations. PhD thesis, Radboud Universiteit Nijmegen, 2005. → pages 4, 12,16, 18[31] G. Frehse. PHAVer: algorithmic verification of hybrid systems past HyTech.International Journal on Software Tools for Technology Transfer, 10(3):263–279, 2008. ISSN 1433-2779. doi:10.1007/s10009-007-0062-x. URLhttp://dx.doi.org/10.1007/s10009-007-0062-x. → pages 4, 12, 16, 18[32] G. Frehse, B. Krogh, and R. Rutenbar. Verifying analog oscillator circuitsusing forward/backward abstraction refinement. In Design, Automation andTest in Europe, 2006. DATE ’06. Proceedings, volume 1, pages 6 pp.–,March 2006. doi:10.1109/DATE.2006.244113. → pages 12[33] G. Frehse, C. L. Guernic, A. Donz, S. Cotton, R. Ray, O. Lebeltel,R. Ripado, A. Girard, T. Dang, and O. Maler. SpaceEx: Scalable verificationof hybrid systems. In G. Gopalakrishnan and S. Qadeer, editors, CAV,volume 6806 of Lecture Notes in Computer Science, pages 379–395.Springer, 2011. ISBN 978-3-642-22109-5. → pages 4, 7, 12, 16, 17[34] S. Gao, S. Kong, and E. Clarke. Satisfiability modulo ODEs. In FormalMethods in Computer-Aided Design (FMCAD), 2013, pages 105–112, Oct2013. doi:10.1109/FMCAD.2013.6679398. → pages 4, 16[35] A. Ghosh and R. Vemuri. Formal verification of synthesized analog designs.In Computer Design, 1999. (ICCD ’99) International Conference on, pages40–45, 1999. doi:10.1109/ICCD.1999.808362. → pages 12[36] R. K. Granlund, T. The GNU multiple precision arithmetic library version4.0. 2001. → pages 1675[37] M. Greenstreet and I. Mitchell. Reachability analysis using polygonalprojections. In Hybrid Systems: Computation and Control, volume 1569 ofLecture Notes in Computer Science, pages 103–116. Springer BerlinHeidelberg, 1999. ISBN 978-3-540-65734-7.doi:10.1007/3-540-48983-5 12. → pages 4, 11, 16, 17[38] M. R. Greenstreet and I. Mitchell. Integrating projections. In T. A.Henzinger and S. Sastry, editors, Proceedings of the First InternationalWorkshop on Hybrid Systems: Computation and Control, pages 159–174,Berkeley, California, Apr 1998. → pages 19[39] C. Guernic and A. Girard. Reachability analysis of hybrid systems usingsupport functions. In Computer Aided Verification, volume 5643 of LectureNotes in Computer Science, pages 540–554. Springer Berlin Heidelberg,2009. ISBN 978-3-642-02657-7. doi:10.1007/978-3-642-02658-4 40. →pages 18[40] S. Gupta, B. Krogh, and R. Rutenbar. Towards formal verification of analogdesigns. In Computer Aided Design, 2004. ICCAD-2004. IEEE/ACMInternational Conference on, pages 210–217, Nov 2004.doi:10.1109/ICCAD.2004.1382573. → pages 12[41] W. Hartong, R. Klausen, and L. Hedrich. Formal verification for nonlinearanalog systems: Approaches to model and equivalence checking. InR. Drechsler, editor, Advanced Formal Verification, pages 205–245. SpringerUS, 2004. ISBN 978-1-4020-7721-0. doi:10.1007/1-4020-2530-0 6. →pages 11[42] L. Hedrich and E. Barke. A formal approach to nonlinear analog circuitverification. In Computer-Aided Design, 1995. ICCAD-95. Digest ofTechnical Papers., 1995 IEEE/ACM International Conference on, pages123–127, Nov 1995. doi:10.1109/ICCAD.1995.480002. → pages 10[43] T. Henzinger. The theory of hybrid automata. In Logic in Computer Science,1996. LICS ’96. Proceedings., Eleventh Annual IEEE Symposium on, pages278–292, Jul 1996. doi:10.1109/LICS.1996.561342. → pages 13, 14[44] T. Henzinger, P.-H. Ho, and H. Wong-Toi. HYTECH: the next generation. InReal-Time Systems Symposium, 1995. Proceedings., 16th IEEE, pages56–65, Dec 1995. doi:10.1109/REAL.1995.495196. → pages 4, 12, 15[45] T. Henzinger, B. Horowitz, R. Majumdar, and H. Wong-Toi. BeyondHyTech: Hybrid systems analysis using interval numerical methods. In76N. Lynch and B. Krogh, editors, Hybrid Systems: Computation and Control,volume 1790 of Lecture Notes in Computer Science, pages 130–144.Springer Berlin Heidelberg, 2000. ISBN 978-3-540-67259-3.doi:10.1007/3-540-46430-1 14. → pages[46] T. Henzinger, J. Preussig, and H. Wong-Toi. Some lessons from theHYTECH experience. In Decision and Control, 2001. Proceedings of the40th IEEE Conference on, volume 3, pages 2887–2892 vol.3, 2001.doi:10.1109/.2001.980714. → pages 4, 12, 15[47] T. A. Henzinger. The theory of hybrid automata. 11th Annual Symposium onLogic in Computer Science (LICS), pages 278–292, 1996. → pages 4[48] A. Jesser and L. Hedrich. A symbolic approach for mixed-signal modelchecking. In Design Automation Conference, 2008. ASPDAC 2008. Asia andSouth Pacific, pages 404–409, March 2008.doi:10.1109/ASPDAC.2008.4483984. → pages 6, 19[49] J. Kim, M. Jeeradit, B. Lim, and M. Horowitz. Leveraging designer’s intent:A path toward simpler analog CAD tools. In Custom Integrated CircuitsConference, 2009. CICC ’09. IEEE, pages 613–620, Sept 2009.doi:10.1109/CICC.2009.5280741. → pages 6[50] K. Kundert. Introduction to rf simulation and its application. Solid-StateCircuits, IEEE Journal of, 34(9):1298–1319, Sep 1999. ISSN 0018-9200.doi:10.1109/4.782091. → pages 6[51] R. Kurshan and K. McMillan. Analysis of digital circuits through symbolicreduction. Computer-Aided Design of Integrated Circuits and Systems, IEEETransactions on, 10(11):1356–1371, Nov 1991. ISSN 0278-0070.doi:10.1109/43.97615. → pages 11[52] A. Kurzhanski and P. Varaiya. Ellipsoidal techniques for reachabilityanalysis: internal approximation. Systems and Control Letters, 41(3):201 –211, 2000. ISSN 0167-6911. doi:10.1016/S0167-6911(00)00059-1. →pages 4[53] T. Larrabee. Test pattern generation using boolean satisfiability.Computer-Aided Design of Integrated Circuits and Systems, IEEETransactions on, 11(1):4–15, Jan 1992. ISSN 0278-0070.doi:10.1109/43.108614. → pages 1077[54] H. Lin, P. Li, and C. Myers. Verification of digitally-intensive analog circuitsvia kernel ridge regression and hybrid reachability analysis. In DesignAutomation Conference (DAC), 2013 50th ACM / EDAC / IEEE, pages 1–6,May 2013. → pages 20[55] J. Lygeros, K. Johansson, S. Simic, J. Zhang, and S. Sastry. Dynamicalproperties of hybrid automata. Automatic Control, IEEE Transactions on, 48(1):2–17, Jan 2003. ISSN 0018-9286. doi:10.1109/TAC.2002.806650. →pages 6[56] L. W. Nagel and D. O. Pederson. Spice (simulation program with integratedcircuit emphasis). Memorandum No. ERL-M382, University of California,Berkeley, Apr. → pages 2[57] S. Osher and J. A. Sethian. Fronts propagating with curvature dependentspeed: Algorithms based on Hamilton-Jacobi Formulations. JOURNAL OFCOMPUTATIONAL PHYSICS, 79(1):12–49, 1988. → pages 4[58] J. Queille and J. Sifakis. Specification and verification of concurrent systemsin CESAR. In M. Dezani-Ciancaglini and U. Montanari, editors,International Symposium on Programming, volume 137 of Lecture Notes inComputer Science, pages 337–351. Springer Berlin Heidelberg, 1982. ISBN978-3-540-11494-9. doi:10.1007/3-540-11494-7 22. → pages 11[59] A. Salem. Semi-formal verification of VHDL-AMS descriptions. In Circuitsand Systems, 2002. ISCAS 2002. IEEE International Symposium on,volume 5, pages 333–336, 2002. doi:10.1109/ISCAS.2002.1010708. →pages 11[60] D. Sauer. Analog to digital converter, nov 1993. US Patent5,262,779,http://www.google.com/patents/US5262779. → pages 2[61] B. Silva, K. Richeson, B. H. Krogh, and A. Chutinan. Modeling andverification of hybrid dynamical system using CheckMate. In ADPM 2000,2000. → pages 4, 12, 16[62] P. Stephan, R. Brayton, and A. Sangiovanni-Vincentelli. Combinational testgeneration using satisfiability. Computer-Aided Design of IntegratedCircuits and Systems, IEEE Transactions on, 15(9):1167–1176, Sep 1996.ISSN 0278-0070. doi:10.1109/43.536723. → pages 10[63] P. Varaiya. Reach set computation using optimal control. In In Proc. KITWorkshop on Hybrid Systems, 1998. → pages 15, 1678[64] Z. Wang, N. Abbasi, R. Narayanan, M. Zaki, G. Al-Sammane, and S. Tahar.Verification of analog and mixed signal designs using online monitoring. InMixed-Signals, Sensors, and Systems Test Workshop, 2009. IMS3TW ’09.IEEE 15th International, pages 1–8, June 2009.doi:10.1109/IMS3TW.2009.5158695. → pages 6, 19[65] J. Wei, Y. Peng, G. Yu, and M. Greenstreet. Verifying global convergencefor a digital phase-locked loop. In Formal Methods in Computer-AidedDesign (FMCAD), 2013, pages 113–120, Oct 2013.doi:10.1109/FMCAD.2013.6679399. → pages iv, 12, 13[66] R. Wilson. Under the lid: Analog test is suddenly the critical ingredient.EDN, 7 January 2010. http://www.edn.com/electronics-news/4312835/Under-the-Lid-Analog-test-is-suddenly-the-critical-ingredient. → pages 3[67] C. Yan. Coho: A verification tool for circuit verification by reachabilityanalysis. Master Dissertation in the University of British Columbia,November, 2006. → pages 4, 12, 16, 17[68] C. Yan. Reachability analysis based circuit-level formal verification. Ph.DDissertation in the University of British Columbia, November 2011. →pages 7, 51[69] C. Yan and M. Greenstreet. Circuit level verification of a high-speed toggle.In Formal Methods in Computer Aided Design, 2007. FMCAD ’07, pages199–206, Nov 2007. doi:10.1109/FAMCAD.2007.39. → pages 12, 17[70] C. Yan and M. Greenstreet. Verifying an arbiter circuit. In Formal Methodsin Computer-Aided Design, 2008. FMCAD ’08, pages 1–9, Nov 2008.doi:10.1109/FMCAD.2008.ECP.11. → pages 12, 17[71] C. Yan, M. Greenstreet, and S. Yang. Verifying global start-up for a Mo¨biusring-oscillator. Formal Methods in System Design, pages 1–27, 2013. ISSN0925-9856. doi:10.1007/s10703-013-0204-6. → pages 12, 17[72] G. Yu. Spectre simulation of a ring-oscillator in a 65nm CMOS process.personal communication, Aug. 2013. → pages 2779Appendix APFD Glitch ProofA.1 ProofTheorem 1. The PFD glitch will not happen if v is in bound (V blow,V bhi). Let :V blow =g2 fref (cmin−ccenter)K + fref cmin− cminKδ ,V bhi =g2 fref (cmax−ccenter)K + fref cmax + cmaxKδ ,θ∆0Because of Reset Delay , signal maybe missed, which cause the phase headingtowards the wrong direction. For this to be happened, the reset delay should belong enough to cover the next uprising edge of the reference signal. This proofgives the condition for the reset and proves that for a given set of parameters, resetwill not happen for some voltage bounds.If fDCO is way smaller than fref , the trajectory will hit the left side in v× cspace(c saturated in mode 4), then keep climbing until θ∆ >= 0 . What we worriedabout the previous version of proof is that, in mode 4, the reset will probably hap-pen, in which case the trajectory will leave mode 4(go right) and oscillate. For thisto be happened, the θ∆ must decrease to some threshold −δ (δ > 0). That is:80Definition.Reset : θ∆0if θ∆ <−δ thenθ∆ = θ∆+2piend ifLemma 1. If the reset happens, then, the following equation must held:−v0− fref × cmincmin×K+g2× fref × (cmin− ccenter)cmin×K2> δProof. Let a = g2× fref (cmin−ccenter)cmin , b =v0cmin− fref , the above equation can be rewroteas:−bK+aK2> δFor the parameter sets Table A.1 that we are using:a =g2 fref (cmin− ccenter)cmin=−0.0022(0.9−1)0.9≈ 0.000444b =v0cmin− fref ,b≤ 0Note that whatever the parameter set we choose in Table A.1, a is a small positivevalue. Now let’s look at the change of θ∆:dθ∆dt= ( fDCO− fref )−Kθ∆= g(t)−Kθ∆(A.1)In mode 4, c is saturated (c = cmin) and v is a linear function of t.g(t) = fDCO− fref =vc− fref =vcmin− fref=v0 +g2 fref (cmin− ccenter)tcmin− fref=(g2 fref (cmin− ccenter)cmin)t +v0cmin− fref= at +b(A.2)81Now combine with Equation A.1 and Equation A.2, we have:dθ∆dt= g(t)−Kθ∆ = at +b−Kθ∆ (A.3)we can solve the above equation and get:θ∆ = (θ∆01−bK+aK2)e−Kt +aKt +(bK−aK2) (A.4)For the above Equation A.4, we are interested in its minimum value. If θ∆0 >−δ , Min(θ∆) <−δ , then reset will happen. Let s =− bK +aK2 , we have s > 0 andθ∆ = (θ∆0 + s)e−Kt − s+aKt (A.5)• As t increase, the linear term aK t will dominate and so θ∆ will be positive.And this is the time when the trajectory leave the left wall (or mode 4 inSPACEEXmodel).• When t is small, aK t is small (a is small and positive, K ∈ (0,1)). θ∆0 >−δ ,– choose θ∆0, so that θ∆0 + s < 0, then,(θ∆0 + s)e−Kt − s+aKt > θ∆0 >−δccenterso, in this case, the reset will not happen.– choose θ∆0, so that θ∆0 + s≥ 0, then,(θ∆0 + s)e−Kt − s+aKt >−s+Min((θ∆0 + s)e−Kt +aKt) >−swhen θ∆0 =−s, above equation has the minimum value −s. But if s isreally big (v0 is big while K is small) and−s <−δ < θ∆0 < 0 then, thereset will happen. And this is equivalent to: s > δ , which is what wewant to prove after substituting s with − bK +aK2 .1initial value of θ∆82So far we have proved the condition of reset. According to Lemma 1, it is easyto see that :if v0 ≥g2 fref (cmin−ccenter)K + fref cmin− cminKδ thenreset will not happenend ifSimilarly, we can analyze the high voltage bound when fDCO is way larger thanthe reference frequency.if v0 ≤g2 fref (cmax−ccenter)K + fref cmax + cmaxKδ thenreset will not happenend ifLet :V blow =g2 fref (cmin−ccenter)K + fref cmin− cminKδ ,V bhi =g2 fref (cmax−ccenter)K + fref cmax + cmaxKδ ,we now prove (by SPACEEX) that for the given set of parameter, V bhi ≥ v ≥V blow will hold after some time t. Which means reset will finally disappear aftervoltage reach the bound [V blow,V bhi].83The current model for the DPLL is:satc = ((c = cmin)∧ (θδ < 0))∨((c = cmax)∧ (θδ > 0))satv = ((v = vmin)∧ (c > ccenter))∨((v = vmax)∧ (c < ccenter))Resetθ∆ = ((θ∆ ≤−2pi)∧ (θ∆ ≥ 2pi)c˙ = −sign(θ∆)g1 fref , if ¬satc= 0, if satcv˙ = g2 fref (c− ccenter), if ¬satv= 0, if satvθ˙∆ = fDCO(c,v)− fref −Kθ∆ if ¬Resetθ∆reset θ∆ to θ∆−2pisign(θ∆), if Resetθ∆(3.6)The parameters used in the experiment is:Table A.1: Four Sets of Parametersparameters Set 1. Set 2. Set 3. Set 4.c [0.9, 1.1] [0.9, 1.1] [0.9, 1.1] [0.9, 1.1]v [1.0, 2.5] [0.2, 3.8] [0.1, 1.5] [1.5, 2.5]δ 5.0 5.0 5.0 5.0fref 2.0 3.0 1.0 2.0ccenter 1.0 1.0 1.0 1.0g1 -0.01 -0.003125 -0.003125 -0.01g2 g1/5 g1/5 g1/5 g1/5K 0.1 0.8 0.8 0.7For parameter set 2 and set 3 the voltage bound is [−0.9,7.7] and [−2.7,5.5]respectively. For parameter set 3 the voltage bound is [−1.3460,6.0460]. Com-pared with their interval of v, there is no need to give further proof for the reset willnever happen. For the parameter set 1, V blow = 1.3540, V bhi = 2.7460. Obviously,we don’t need to worry about the high bound because of the interval of v given84above. Therefore, next, we show that v will increase to V blowAs we can see from equation (3.6), if c < ccenter, v will keep increase. So wejust need to make sure that when reset happened and v∈ [v0,V blow], the interval forc is less than ccenter. To avoid oscillation in SPACEEX, which will probably causeexplosion, we make v to be an interval and divide them into 7 regions, this will alsofacilitate the calculation of fDCO. The SPACEEXmodel is:c˙ = sign(θ∆)−g2 frefθ˙∆ = fDCO− fref −Kθ∆∈ p1interval(v)− p2c+ p3− fref −Kθ∆⊕〈E3 fref 〉(A.6)The largest value of c depends on the time when θ∆ hit 0 after reset. For ourexperiment, after reset, θ∆ = −5+ 2 ∗ pi ≈ 1.2832, so t is the time for θ∆ to dropfrom 1.2832 to 0. However, in real circuit, the detection of θ∆ to 0 could happenedanywhere to the next cycle. So we also model this into SPACEEX by introducingan interval for θ∆ during transition.Figure A.1: Change of c in Region One, Four and SevenFigure A.1 shows the change of c overtime, value of c is clearly less than ccenter,indicating that v will increase to V b and reset will disappear finally. Also, for asthe value of v increase( region number increase from one to seven), the value of calso increased. Because, the large the v, the smaller the difference of | fDCO− fref |,85which would cause more time for θ∆ to decrease from (2pi− δ ) to 0 and cause cto go even further from cmin.86Appendix BModel Derivation for fDCO and Z3ProofThis appendix derives a COHO-REACHmodel with variable shifting such that vand c are centered at the line v/c = fref in v/c space. An error term is included asdifferential an inclusion to produce a COHO-REACHmodel. We use the Z3 SMTsolver [25] to verify this derivation.First, recall that the dynamics of the digital PLL are described as a set of dif-ferential equations as derived in Section 3.2.2:satc = ((c = cmin)∧ (θδ < 0))∨((c = cmax)∧ (θδ > 0))satv = ((v = vmin)∧ (c > ccenter))∨((v = vmax)∧ (c < ccenter))Resetθ∆ = ((θ∆ ≤−2pi)∧ (θ∆ ≥ 2pi)c˙ = −sign(θ∆)g1 fref , if ¬satc= 0, if satcv˙ = g2 fref (c− ccenter), if ¬satv= 0, if satvθ˙∆ = fDCO(c,v)− fref −Kθ∆ if ¬Resetθ∆reset θ∆ to θ∆−2pisign(θ∆), if Resetθ∆(3.6)87Figure B.1: Variable Shifting!!,!"!!!!!,!!!!!!!,!"!!! !!,!!!!!!!,!!!!!!,!!!!(!!,!!!!!,!!!)!∆!! !∆!! !Where,fDCO(c,v) = vc (3.5)Let [vi,lo,vi,hi] be the voltage bounds for the ith such rectangle and [ci,lo,ci,hi] bethe capacitance bounds. For convenience, I’ll define:vi,0 = 12(vi,hi + vi,lo) ∆vi =12(vi,hi− vi,lo)ci,0 = 12(ci,hi + ci,lo) ∆ci =12(ci,hi− ci,lo)fi,0 =vi,0ci,0(B.1)Consider the case where v ∈ [vi,lo,vi,hi] and c ∈ [ci,lo,ci,hi].v1,i = v− vi,0 c1,i = c− ci,0 (B.2)Figure B.1 shows the above variable shifting. To simplify the derivations, it’sconvenient to normalize voltage, capacitance, and frequency to their midpoint val-ues. In particular, let:v˜i =v1,ivi,0∆v˜i =∆vivi,0c˜i =c1,ivi,0∆c˜i =∆cici,0f˜i =1+ v˜i1+ c˜i(B.3)88For the digital PLL, the nonlinear term mainly comes from fDCO. A linearinclusion for fDCO can be readily obtained given an inclusion for f˜i. We have:fDCO =vc, Eq. 3.5= fi,0 f˜i, Eqs. B.2 and B.3f˜i = (1+ v˜i− c˜i− v˜ic˜i)(1+c˜2i1− c˜2i)(B.4)The non-linearities are in the v˜ic˜i term of the first parenthesized factor, and thec˜2i /(1− c˜2i ) term of the second. Thus, the errors diminish quadratically with thelengths of the sides of the rectangles used to partition the v× c space, as we wouldhope. We now bound each of these terms.Clearly,v˜ic˜i ∈ 〈∆v˜i∆c˜i〉, (B.5)c˜2i1− c˜2i∈[0,∆c˜2i1−∆c˜2i](B.6)LetEi,1 = ∆v˜i∆c˜iEi,2 =∆c˜2i2(1−∆c˜2i )Ai = 1+Ei,2Ei,3 = AiEi,1 +Ei,2 +(∆v˜i +∆v˜i∆c˜i +∆c˜i)Ei,2= Ei,1 +(1+∆v˜i +∆c˜i +2Ei,1)Ei,2Erri = Ei,3 fi,0(B.7)For v ∈ [vi,lo,vi,hi] and c ∈ [ci,lo,chii], we have1+ v˜i− c˜i− v˜ic˜i ∈ 1+ v˜i− c˜i⊕〈Ei,1〉1+ c˜2i1−c˜2i∈ Ai⊕〈Ei,2〉f˜i ∈ (1+ v˜i− c˜i⊕〈Ei,1〉)(Ai⊕〈Ei,2〉)⊆ Ai(1+ v˜i− c˜i)⊕〈Ei,3〉fDCO ∈ Ai(1+ v˜i− c˜i) fi,0⊕〈Ei,3 fi,0〉(B.8)89Let c1 = c− c0, v1 = v− v0. The linearized differential model can be writtenas:c˙1 =−sign(θ∆)g1 frefv˙1 ∈ g2 fref c1−g2 fref (ccenter−vmidfref)⊕g2vmid v˜θ˙∆ ∈ A(1+v1vmid−c1cmid)fref − fref −Kθ∆⊕〈E3 fref 〉A = 1+E2E1 = c˜v˜E2 =c˜22(1− c˜2)E3 = E1 +(1+ v˜+ c˜+2E1)E2c˜ =chi− clochi + clov˜ =vhi− vlovhi + vlo(B.9)B.1 Z3 Validation of the Derived ModelThe verification of global convergence for the digital PLL depends on the correct-ness of the linearization presented above. This was a manual derivation. Thus, wesought a way to prove the correctness of the results. To do so, we used the Z3 SMTsolver. Initially, we attempted simply wrote the final linear inclusion for fDCO, andasked Z3 to prove that it contains the non-linear formula. Z3 failed. With someexperiments, we concluded that the many definitions and constraints caused Z3 totake an unfruitful path when subdividing the space while trying to verify the claim.To avoid this problem, we wrote a simple theorem prover on top of Z3s pythonAPI. Basically, this allows us to prove lemmas from the original constraints or pre-viously proven lemmas. For each such proof, we state what lemmas or hypothesesto use. This allows us to make sure that Z3 is given the relevant constraints foreach lemma. With this approach, verifying linear inclusion for fDCO was relatively90straightforward. The python code is presented in the remainder of this appendix.91from z3 impor t ∗c lass TheoremProver :c lass Proo fFa i l u re ( Except ion ) :def i n i t ( se l f , c laim , why ) :s e l f . c la im = cla ims e l f . why = why# end i n i tdef ve rbos i t y ( s e l f ) : r e t u rn 1def s t r ( s e l f ) :s = ” Fa i led to es t ab l i s h clause ”i f ( s e l f . ve rbos i t y ( ) > 0 ) :s = s + ” :\ n clause = ” + s e l f . c la im . s t r ( )s = s + ”\n because ” + s e l f . whyre tu rn s# end s t r# end Learn ingFa i lu redef i n i t ( se l f , hypotheses ={} ) :s e l f . lemmas = hypothesess e l f . f a i l u r e s = {}def prove ( se l f , c laim , name, hypotheses = [ ] ) :i f ( not (name i s None ) and s e l f . lemmas . c on t a i n s (name ) ) :msg = ” a lemma / hypothes is w i th name ” + repr (name) + ” a l ready ex i s t s ”r a i se TheoremProver . P roo fFa i l u re ( claim , msg)s = Solver ( )f o r h i n hypotheses [ : ] :i f ( s e l f . lemmas . c on t a i n s ( h ) ) :s . add ( s e l f . lemmas [ h ] )e lse :r a i se TheoremProver . P roo fFa i l u re ( claim , ” unknown lemma ” + repr ( h ) )s . add ( Not ( c la im ) )s ta tus = s . check ( )i f ( s ta tus == unsat ) :i f ( not (name i s None ) ) :s e l f . lemmas [name ] = c la imre tu rn Truee l i f ( s ta tus == sat ) :s e l f . f a i l u r e s [ name ] = ( hypotheses , claim , s . model ( ) )p r i n t ” counterexample ”p r i n t s . model ( )r e t u rn Falsee lse :92r e t u rn s ta tusdef cex ( se l f , name ) :r e t u rn s e l f . f a i l u r e s [ name ]from z3 impor t ∗from TheoremProver impor t TheoremProverc lass ProofFa i ledExcept ion ( Except ion ) :def i n i t ( se l f , why ) :s e l f . why = whydef s t r ( s e l f ) :r e t u rn s e l f . whydef hooray (name ) :p r i n t ” Hooray −− we proved ” , namedef myProve ( prover , claim , name, hypotheses ) :s ta tus = prover . prove ( claim , name, hypotheses )i f ( s ta tus == True ) :hooray (name)re t u rn Truee l i f ( s ta tus == False ) :r a i se ProofFa i ledExcept ion (name + ” : re fu ted ( counterexample generated ) ” )e lse :r a i se ProofFa i ledExcept ion (name + ” : i nconc lus i ve ” )def prove er r bnds ( ) :# dec lare va r i ab l es f o r ( normal ized ) capaci tance[ c , c lo , c h i ] = Reals ( [ ” c ” , ” c l o ” , ” c h i ” ] )[ c 0 , dc , c 1 , cc , dcc ] = Reals ( [ ” c 0 ” , ” dc ” , ” c 1 ” , ” cc ” , ” dcc ” ] )c hyps = { ”0 < c l o ” : 0 < c lo ,” c . bnds ” : And ( c l o <= c , c <= c h i ) ,” c 0 . def ” : c 0 == ( c h i + c l o ) / 2 ,” dc . def ” : dc == ( c h i − c l o ) / 2 ,” c 1 . def ” : c 1 == c − c 0 ,” cc . def ” : cc == c 1 / c 0 ,” dcc . def ” : dcc == dc / c 0 }# dec lare va r i ab l es f o r ( normal ized ) capaci tance[ v , v lo , v h i ] = Reals ( [ ” v ” , ” v l o ” , ” v h i ” ] )[ v 0 , dv , v 1 , vv , dvv ] = Reals ( [ ” v 0 ” , ” dv ” , ” v 1 ” , ” vv ” , ” dvv ” ] )v hyps = { ”0 < v l o ” : 0 < v lo ,” v . bnds ” : And ( v l o <= v , v <= v h i ) ,” v 0 . def ” : v 0 == ( v h i + v l o ) / 2 ,93” dv . def ” : dv == ( v h i − v l o ) / 2 ,” v 1 . def ” : v 1 == v − v 0 ,” vv . def ” : vv == v 1 / v 0 ,” dvv . def ” : dvv == dv / v 0 }# dec lare va r i ab l es f o r ( normal ized ) frequency[ f DCO , f 0 , ff DCO , f5 ] = Reals ( [ ” f DCO ” , ” f 0 ” , ” ff DCO ” , ” f5 ” ] )f hyps = { ” f DCO . def ” : f DCO == v / c ,” f 0 . def ” : f 0 == v 0 / c 0 ,” ff DCO . def ” : ff DCO == f DCO / f 0 ,” f5 . def ” : f5 == (1 + vv − cc − vv∗cc )∗ (1 + ( cc∗cc ) / ( 1 − ( cc∗cc ) ) ) }# some quan t i t i e s used in c a l c u l a t i n g e r r o r bounds[ cx , dcx , A, e1 , e2 , e3 , errbnd ] = Reals ( [ ” cx ” , ” dcx ” , ”A” , ” e1 ” , ” e2 ” , ” e3 ” , ” errbnd ” ] )x hyps = { ” cx . def ” : cx == ( cc∗cc ) / (1− ( cc∗cc ) ) ,” dcx . def ” : dcx == ( dcc∗dcc ) / (1− ( dcc∗dcc ) ) ,”A . def ” : A == 1 + e2 ,” e1 . def ” : e1 == dvv∗dcc ,” e2 . def ” : e2 == dcx /2 ,” e3 . def ” : e3 == e1 + (1 + dvv + dcc + 2∗e1)∗e2 ,” errbnd . def ” : errbnd == e3∗ f 0 }# create the theorem proverhyps = d i c t ( d i c t ( d i c t ( d i c t (∗∗ c hyps ) , ∗∗v hyps ) , ∗∗ f hyps ) , ∗∗x hyps )p = TheoremProver ( hyps )t r y :# bounds on c 0 , dcc , v 0 , dvv , f 0myProve (p , 0 < c 0 , ”0 < c 0 ” , c hyps . keys ( ) )myProve (p , And(0 <= dcc , dcc < 1) , ” dcc . bnds ” , c hyps . keys ( ) )myProve (p , And(−dcc <= cc , cc <= dcc ) , ” cc . bnds ” , c hyps . keys ( ) )myProve (p , 0 < 1 − cc∗cc , ”0 < 1 − cc∗cc ” , [ ” dcc . bnds ” , ” cc . bnds ” ] )myProve (p , 0 < 1 − dcc∗dcc , ”0 < 1 − dcc∗dcc ” , [ ” dcc . bnds ” ] )myProve (p , 0 < v 0 , ”0 < v 0 ” , v hyps . keys ( ) )myProve (p , And(0 <= dvv , dvv < 1) , ” dvv . bnds ” , v hyps . keys ( ) )myProve (p , And(−dvv <= vv , vv <= dvv ) , ” vv . bnds ” , v hyps . keys ( ) )myProve (p , 0 < f 0 , ”0 < f 0 ” , [ ” f 0 . def ” , ”0 < c 0 ” , ”0 < v 0 ” ] )# c o l l e c t commonly used combinat ions o f lemmas toge therpos lems = [ ” 0 < c l o ” , ”0 < c 0 ” , ”0 < v l o ” , ”0 < v 0 ” , ”0 < v 0 ” ,”0 < 1 − cc∗cc ” , ”0 < 1 − dcc∗dcc ” , ”0 < f 0 ” ]cdef lems = [ ” c 0 . def ” , ” c 1 . def ” , ” cc . def ” ]vdef lems = [ ” v 0 . def ” , ” v 1 . def ” , ” vv . def ” ]fde f lems = [ ” f DCO . def ” , ” f 0 . def ” , ” ff DCO . def ” , ” f5 . def ” ]def lems = cdef lems + vdef lems + fde f lems94# check equat ion 5myProve (p , ff DCO == f5 , ” eq5 ” , pos lems + def lems + c hyps . keys ( ) + v hyps . keys ( ) )myProve (p , f DCO == f5∗ f 0 , ” eq5b ” , pos lems + [ ” eq5 ” ] + f hyps . keys ( ) )# check equat ion 7myProve (p , And(0 <= cx , cx <= dcx ) , ” eq7 ” , pos lems + c hyps . keys ( ) + x hyps . keys ( ) )# check equat ion 9c c s t u f f = [ ” dcc . bnds ” , ” cc . bnds ” , ”0 < 1 − cc∗cc ” , ”0 < 1 − dcc∗dcc ” ]v v s t u f f = [ ” dvv . bnds ” , ” vv . bnds ” ]# not sure why we need eq9a to prove eq9bmyProve (p , And (A∗(1 + vv − cc ) − e3 <= ff DCO , ff DCO <= A∗(1 + vv − cc ) + e3 ) ,” eq9a ” , pos lems + [ ” eq5 ” , ” eq7 ” ] + c hyps . keys ( ) + v hyps . keys ( ) + f hyps . keys ( )+ x hyps . keys ( ) )myProve (p , And (A∗(1 + vv − cc ) − e3 <= f5 , f5 <= A∗(1 + vv − cc ) + e3 ) , ” eq9b ” , pos lems+ [ ” eq7 ” , ” eq9a ” ] + c c s t u f f + v v s t u f f + fde f lems + x hyps . keys ( ) )myProve (p , And (A∗(1 + vv − cc )∗ f 0 − e3∗ f 0 <= f DCO , f DCO <= A∗(1 + vv − cc )∗ f 0 + e3∗ f 0 ) ,” eq9 ” , [ ” eq5b ” , ” eq9a ” ] + pos lems + fde f lems )p r i n t ” A l l c la ims proven ! ”except ( ProofFa i ledExcept ion , TheoremProver . P roo fFa i l u re ) as oops :impor t sysp r i n t >> sys . s tde r r , oopsre t u rn p95Appendix CModel Derivation for SpaceexModel 2This appendix derives the spaceex model 2. In model 2, fDCO is close to fref , θ∆frequently changes its sign, causing a large number of mode transitions that wouldobscure the progress of v in the SPACEEX analysis. To show that the value of fDCOkeeps making progress, we use a new parameter w, where:w = fDCO(c,v)− frefw˙ = ddt fDCO(c,v)− f˙refAssume that f˙ref = 0 and concentrate on the fDCO term:fDCO(c,v) = vc , Equation 3.5fDCO = v˙c −vc2 c˙, calculusc˙ = −sign(θ∆)g1 fref , Equation 3.6v˙ = g2 fref (c− ccenter), Equation 3.6It would be sound to use the inclusions for c˙ and v˙, but the c and c2 terms in thedenominator of the equation for fDCO must still be handled. Now we have:96fDCO = v˙c −vc2 c˙, shown above= g2 fref (c−ccenter)c +sign(θ∆)g1 fref vc2 , Equation 3.6= fref(g2−g2ccenterc +sign(θ∆)g1vc2), algebra= fref(g2−g2ccenterc +sign(θ∆)g1 fDCOc), Equation 3.5= fref(g2−g2ccenterc +sign(θ∆)g1( fref +w)c), w = fDCO− fref= fref(g2 +g1sign(θ∆) fref−g2ccenterc +sign(θ∆)g1c w), algebraNext, get the inclusion that we need:fref g2 is a constant.fref (g1sign(θ∆) fref −g2ccenter)/c can be bounded by an interval depending on thesign(θ∆). For all parameter sets from Table 4.1, g1 fref < g2ccenter < 0, c > 0,fref > 0, and g1 < 0. If θ∆ > 0, then fref (g1sign(θ∆) fref − g2ccenter)/c < 0,and the bounding interval for this term is fref (g1sign(θ∆) fref −g2ccenter)[1/clo,1/chi].If θ∆≤ 0, then the bounding interval is fref (g1sign(θ∆) fref −g2ccenter)[1/chi,1/clo].For both of these, clo and chi are bound on c for the mode. For example, inModel 2lo, we could let clo = cmin and chi = ccenter. If we subdivide Mode2lo into smaller intervals, then we can choose clo and chi accordingly. Let I1denote this interval. Writing out the equations:I1 = (g1sign(θ∆) fref −g2ccenter)[1/clo,1/chi], if θ∆ > 0= (g1sign(θ∆) fref −g2ccenter)[1/chi,1/clo], if θ∆ ≤ 0g1sign(θ∆) fref−g2ccenterc ∈ I1fref sign(θ∆)g1c w is nearly the linear term in w that we would like. For the term c inthe denominator. Letchmean =2chiclochi + clo, harmonic meanηc =chi− clo2chiclo97Then,fref sign(θ∆)g1c w ∈fref sign(θ∆)g1chmeanw⊕〈 fref g1ηcwmax〉where wmax is an upper bound on |w|. LetI2 = 〈g1ηcwmax〉and concludesign(θ∆)g1c w ∈sign(θ∆)g1chmeanw⊕ I2Therefore, we get:w˙ ∈ fref(g2 +sign(θ∆)g1chmeanw⊕ I1⊕ I2)98Appendix DCoho Verification CodeThe MATLAB code for COHO-REACH verification is presented as follows:% th i s i s the top f unc t i on to v e r i f y PLL , i t cons i s t s o f f ou r steps ,% which w i l l take the whole v a l i d space and move space i n t o a% narrow s t r i p according to the dynamic o f the p l l% step1 : any i n i t i a l cond i t i on leads to to c=c lo , c=c h i , or the ta=0% step2 : t r a j e c t o r i e s on the c−wal l s even tua l l y leave% step3 : successive cross ings o f the ta=0 move occur w i th f dco% converging to N∗ f r e f% step4 : z ig−zags i n a narrow s t r i p near f dco=N∗ f r e f converge to% c near c midf unc t i on isConverge = P l lV ( p , ds )i f ( nargin<1)p = getP ( ) ;e lsep = getP ( p ) ; % a func t i on t ha t can be used to get a l l parametersendi f ( nargin<2)ds = 0.03;% ds def ines the width o f the narrow s t r i pendisConverge = 1;%step one : t r a j e c t o r i e s s t a r t i n g anywhere i n ( c , v , the ta ) even tua l l y% h i t a wa l l ( i . e c = c min , c = c max ) or the ta crosses 0 . We don ’ t% check f o r v h i t t i n g a wa l l −− should we? No, l e t ’ s not care about i t% f i r s t% Because we a l low a r b i t r a r y i n i t i a l cond i t ions , c can h i t a wa l l w i th% a r b i t r a r y values f o r v and the ta . L ikewise the ta can cross 0 wi th% a r b i t r a r y values f o r c and v . So , we don ’ t need to record the end% regions f o r t h i s step , we j u s t need to know tha t one of these99% hyperplanes i s even tua l l y crossed .data1=step one ( p ) ;i f ( data1 . pass )d isp ( ’ pass step one ’ ) ;e lseisConverge = 0;d isp ( ’ f a i l i n step one ’ ) ;end%step two : t r a j e c t o r i e s on the c=c min or c=c max wa l l even tua l l y% leave the wa l l ( i . e . f DCO crosses N∗ f r e f ) .% The p recond i t i on f o r s tep two i s c=c min ( or c=c max ) and% any values f o r v and the ta .% The post−cond i t i on i s the ta=0 and ( d / d t ) the ta > 0 ( f o r the c=c min% case ) . Do we need to es t ab l i s h anyth ing about v?% we need to show tha t v w i l l even tua l l y en ter the narrow s t r i pdata2=step two (p , ds ) ;the ta2 = abs ( data2 . the ta ( 1 ) ) ;i f ( theta2<abs ( data2 . the ta ( 2 ) ) )the ta2 = abs ( data2 . the ta ( 2 ) ) ;endi f ( data2 . pass )d isp ( ’ pass step two , and wi th the f o l l ow i ng the ta bound ’ ) ;d isp ( data2 . the taL ) ;d isp ( data2 . thetaH ) ;e lseisConverge = 0;d isp ( ’ f a i l i n step two ’ ) ;end% step three : t r a j e c t o r i e s t ha t cross the ta=0 wi th f DCO < N∗ f r e f% w i l l have t h e i r next c ross ing o f the ta=0 ( i n the other d i r e c t i o n )% wi th f DCO c lose r to N∗ f r e f . This DPLL has f a i r l y l a rge time−gain% so , f DCO converges to a value c lose to N∗ f r e f qu i t e qu i c k l y .% As wi th step 2 , we need to show bounds on c and v (when the ta crosses% 0) t ha t s a t i s f y the pre−cond i t i ons f o r step fou r −− step fou r w i l l be% added ” soon ” .data3 1 = s tep th ree (0 , ds , p ) ;data3 2 = s tep th ree (1 , ds , p ) ;i f ( data3 1 . pass && data3 2 . pass )d isp ( ’ pass step three ’ ) ;e lseisConverge = 0;d isp ( ’ f a i l i n step three ’ ) ;end100data4 = s tep fou r ( p , ds , [ 0 , the ta2 ] ) ;i f ( data4 . pass )d isp ( ’ pass step four ’ ) ;e lseisConverge = 0;d isp ( ’ f a i l i n step four ’ ) ;end% isConverge means a l l th ree steps passed . This means t ha t% f DCO i s ” c lose ” to N∗ f r e f . From there , t r a j e c t o r i e s w i l l% zig−zag to b r ing c c lose to cc . We don ’ t check the zig−zags% here . I have another .m f i l e f o r tha t , and I ’ l l add those checks% to t h i s f unc t i on l a t e r .i f ( isConverge ==0)d isp ( ’ the PLL f a i l to converge ! ’ )e lsed isp ( ’ the PLL converges ! ’ )endend % Pl lV% t h i s f unc t i on i s used to check t ha t s t a r t i n g anywhere i n the% va l i d space , i t w i l l e i t h e r end up wi th h i t t i n g the wa l l or the ta = 0func t i on data=step one ( p )data . pass = 1;i f ( nargin<1)p = getP ( ) ; %get a l l the parameterse lsep = getP ( p ) ;end% set de f au l t s f o r PLL model parametersi n i t b b o x = c e l l ( 1 , 0 ) ; k = 1 ;% lower pa r t i n mode1i n i t b b o x{k} = [ p . cmin p . cmax ; . . .p . vmin p . f r e f ∗(p . cmin+p . cmax ) / 2 ; 0 2∗p i ] ; k=k+1;% lower pa r t i n mode2i n i t b b o x{k} = [ p . cmin p . cmax ; . . .p . vmin p . f r e f ∗(p . cmin+p . cmax)/2;−2∗ p i 0 ] ; k=k+1;% h ige r pa r t i n mode1i n i t b b o x{k} = [ p . cmin p . cmax ; . . .p . f r e f ∗(p . cmin+p . cmax ) / 2 p . vmax ;0 2∗p i ] ; k=k+1;% higher pa r t i n mode2i n i t b b o x{k} = [ p . cmin p . cmax ; . . .101p . f r e f ∗(p . cmin+p . cmax ) / 2 p . vmax;−2∗p i 0 ] ;dim = 3; planes = [ 1 , 2 ; 1 , 3 ; 2 , 3 ] ;f o r i = 1 : ki n i t P h = ph createByBox ( dim , planes , i n i t b b o x{ i } ) ;i n i t P h = ph conver t ( i n i tPh , ’ convex ’ ) ;i f ( i /2==1)data1 = Run mode1 ( i n i tPh , [ ] , p ) ;% data would be 1 i f i t s topse lsedata1 = Run mode2 ( i n i tPh , [ ] , p ) ;endi f ( data1 . stop ==0)%we are i n t e r es t ed i n whether i t s tops or notdata . pass = 0;d isp ( ’ f a i l i n step one ! ’ )endendi f ( data . pass ˜=0)data . pass = 1;d isp ( ’ pass step one ! ’ ) ;endend% t h i s f unc t i on i s used to check% 1. f o r t r a j e c t o r i e s i n the wa l l ( l e f t , or r i g h t ) ,% i t w i l l f i n a l l y leave the wa l l% 2. the value o f v when i t leave the wa l l to make sure% tha t i t w i l l en ter the narrow s t r i p def ined by ds%func t i on data=step two (p , ds )i f ( nargin<1)p = getP ( ) ;e lsep = getP ( p ) ;endi f ( nargin<2)ds=0.03;enddata . pass = 1;% f i r s t run mode4 : the ta i n [0 , 2∗p i ] , c = c max% v i s a b i t t r i c k y −− I wanted to assume tha t v / c max > N∗ f r e f . That102% corresponds to the case where v c l imbs down the wa l l to get to the% desi red frequency . Then , my crazy superv iso r po in ted out i n i t i a l% cond i t i ons where theta>0, c ( 0 ) i s j u s t s l i g h t l y smal le r than c max , and% v (0 ) < N∗ f r e f ∗c max . Such a t r a j e c t o r y w i l l h i t the c=c max wa l l . v% w i l l drop ( because c > cc ) . The t r a j e c t o r y w i l l go back to the i ns i de% of the c−v rec tang le . Because v / c < N∗ f r e f , c w i l l decrease and we ’ l l% approach lock from below . I need to handle t h i s case or I ’ l l never% graduate −− ( not r ea l l y , my superv iso r i s n i ce r than t ha t ) .h i v = p . f r e f ∗(p . cmax + ds ) ; % high v bound f o r the narrow s t r i plow v = p . f r e f ∗(p . cmin − ds ) ; % low v bound f o r the narrow s t r i p%def ined the i n i t i a l reg iondim = 2; planes = [ 1 , 2 ] ;% set de f au l t s f o r PLL model parametersi n i t b b o x = c e l l ( 1 , 0 ) ; k = 1 ;i n i t b b o x{k} = [ p . vmin p . cmin∗p . f r e f ;−2∗p i 0 . 0 ] ; k = k+1;% run mode3 , check vi n i t b b o x{k} = [ p . cmin∗p . f r e f , p . vmax;−2∗pi , 0 . 0 ] ; k = k+1;% run mode3i n i t b b o x{k} = [ p . cmax∗p . f r e f p . vmax ;0 2∗p i ] ; k = k+1;%run mode4 , check vi n i t b b o x{k} = [ p . vmin , p . cmax∗p . f r e f ; 0 . 0 , 2∗p i ] ;%run mode4f o r i = 1 : ki n i t P h = ph createByBox ( dim , planes , i n i t b b o x{ i } ) ;i n i t P h = ph conver t ( i n i tPh , ’ convex ’ ) ;i f ( i<3)data2 = Run mode3 (p , i n i t P h ) ;e lsedata2 = Run mode4 (p , i n i t P h ) ;end%the f o l l ow i ng cond i t i on check t ha t t r a j e c t o r i e s w i l l leave wa l l and% enter the s t r i p ( by checking v bound )i f ( data2 . leavewa l l ==0 | | ( i ==1 && data2 . phin i t mode1 . bbox(2 ,1)< low v ) . . .| | ( i ==3 && data2 . phin i t mode2 . bbox(2 ,2)> h i v ) )data . pass = 0;endendi f ( data . pass ˜=0)d isp ( ’ pass i n step two ’ )e lsed isp ( ’ f a i l i n step two ’ )103endend% the aim of t h i s f unc t i on i s used to v e r i f y t ha t :% s t a r t i n g from anywhere i n the va l i d space wi th theta<0% i t w i l l leave the s t r i p and move upwards or downwards% u n t i l i t i s c lose enough to the l i n e v / c = f r e f%% inpu t : a , i n d i c a t i n g the lower or the h igher pa r t o f the space% f o r a = 0 , i n i t i a l reg ion i s the lower par t , a = 1 ,% upper pa r t ds , a parameter i n d i c a t i n g the d is tance% to v / c = f r e f% output : data%% note t ha t f o r each of the s t r i p i t w i l l run u n t i l% ph empty and check containmentf unc t i on data= s tep th ree ( a , ds , p )data . pass=1;i f ( nargin<1)a = 0;endi f ( nargin<2)ds = 0 .03 ;endi f ( nargin<3)p = getP ( ) ;e lsep = getP ( p ) ;end% depending on which pa r t you are working on , de f ine a i n i t i a ldim = 3; planes = [ 1 , 2 ; 1 , 3 ; 2 , 3 ] ;bbox = [ ] ;% For the lower l i ne ,% we know v>=(c−ds )∗ f r e f => p . vmin>=(p . cmax−ds )∗p . f r e f% there fo re , ds low = ( p . cmax∗p . f r e f−p . vmin ) / p . f r e f% Likewise , f o r the h igher par t ,% v <= ( c+ds )∗ f r eq=>p . vmax<=(p . cmin+ds )∗p . f r e f% there fo re , ds h i = ( p . vmax−p . cmin∗p . f r e f ) / p . f r e fds low = ( p . cmax∗p . f r e f−p . vmin ) / p . f r e f ;% the f a s t e s t d is tance we can go downwardsds h i = ( p . vmax−p . cmin∗p . f r e f ) / p . f r e f ;% the f a s t e s t d is tance we can go upwardsstr ipNum = 5; % number o f s t r i p s104s t r ipBound low = ( ds .∗ ( str ipNum : ( −1 ) : 0 ) + . . .ds low .∗ ( 0 : ( + 1 ) : str ipNum ) ) . / str ipNum ;s t r i pBound h i = ( ds .∗ ( str ipNum : ( −1 ) : 0 ) + . . .ds h i .∗ ( 0 : ( + 1 ) : str ipNum ) ) . / str ipNum ;f o r k = str ipNum :−1:1 % check each of the s t r i pi f ( a ) % a ==1 , h igher pa r tds1 =s t r i pBound h i ( k+1) ; ds2=s t r i pBound h i ( k ) ; %ds1>ds2%v >= ( c+ds )∗ f r eq => c − v / f r eq <= −dsi n i t l o w = l p c r ea t e ( [1 ,−1/p . f r e f ,0] ,−ds2 ) ; % >=%v <= ( c+ds )∗ f r eq => −c + v / f r eq <= dsi n i t h i = l p c r ea t e ( [−1 ,1/p . f r e f , 0 ] , ds1);%<=i n i t = lp and ( i n i t l ow , i n i t h i ) ;bbox = [ p . cmin , p . cmax ; p . vmin , p . vmax ; 0 , 2∗p i ] ;e lse % a ==0 , lower pa r tds1 =st r ipBound low ( k+1) ; ds2=st r ipBound low ( k );%ds1>ds2%v >= ( c−ds )∗ f r eq => c − v / f r eq <= dsi n i t l o w = l p c r ea t e ( [1 ,−1/p . f r e f , 0 ] , ds1);%>=%v <= ( c−ds )∗ f r eq => −c + v / f r eq <= −dsi n i t h i = l p c r ea t e ( [−1 ,1/p . f r e f ,0] ,−ds2);%<=i n i t = lp and ( i n i t l ow , i n i t h i ) ;bbox = [ p . cmin , p . cmax ; p . vmin , p . vmax ; −2∗pi , 0 ] ;end% create the i n i t i a l ph and i n v a r i a n t by using the% l i n e a r program of the s t r i pi n i t = lp and ( i n i t , lp createByBox ( bbox ) ) ;%i n i t Ph = ph conver t ( ph createByLP ( dim , planes , i n i t ) , ’ convex ’ ) ;i n i t P h = ph createByLP ( dim , planes , i n i t ) ;%ph d isp lay ( i n i tPh , [ ] , [ ] , [ ] , ’ b ’ ) ; drawnow ;data = [ ] ;ca l lBacks . exi tCond = ha ca l lBacks ( ’ exitCond ’ , ’ phempty ’ ) ;%run u n t i l phase emptyca l lBacks . s l iceCond = @( i n f o ) ( 0 ) ;i f ( a ) %a==1ha = mode1 ( [ ] , ca l lBacks , i n i tPh , i n i t , p ) ;da ta f i l e name = [ ha . name, ’ ’ , ha . snames{1} , ’ . mat ’ ] ;reachData = LoadData ( da ta f i l e name ) ;ph0 = reachData{end−1};e lse %a==0ha = mode2 ( [ ] , ca l lBacks , i n i tPh , i n i t , p ) ;da ta f i l e name = [ ha . name, ’ ’ , ha . snames{1} , ’ . mat ’ ] ;reachData = LoadData ( da ta f i l e name ) ;ph0 = reachData{end−1};endres = checkPass ( ph0 , p . f r e f , a ) ;r = [ ds2 , ds1]>=res ;105%the lower bound may over lap , the upper bound of%the s t r i p s t r i c t l y moves outi f ( r ( 2 ) )d isp ( ’ success fu l l y leave s t r i p ’ ) ;e lsed isp ( ’ f a i l to leave s t r i p ’ ) ;data . pass= 0;end%ph d isp lay ( reachData . sets{end− 1} , [ ] , [ ] , [ ] , ’ r ’ ) ; drawnow ;endi f ( data ==0)d isp ( ’ f a i l ’ )e lsedata =1;d isp ( ’ pass ’ )endend% t h i s f unc t i on run step four , where we are i n the narrow s t r i p% we use the c1−v1 model , f o r more d e t a i l s o f the de r i va t i on ,% p ls see my thes is , chapter 5 .f unc t i on data = s t ep f ou r ( p , ds , the ta )i f ( nargin<1)p = getP ( ) ;e lsep = getP ( p ) ;endi f ( nargin<2)ds = 0 .03 ;endi f ( nargin<3)the ta = [0 0 . 5 ] ;enddata . pass = 1;mid lo = p . f r e f ∗(p . cc − ds ) ; % mid boxmid h i = p . f r e f ∗(p . cc + ds ) ;l o l o = p . f r e f ∗(p . cmin − ds);% l e f th i h i = p . f r e f ∗(p . cmax + ds ) ;NumI = 5; %number o f V i n t e r v a l sLoVI = ( l o l o .∗ (NumI : ( −1) :0 ) + mid lo .∗ ( 0 : ( + 1 ) : NumI ) ) . / NumI ;HiVI = ( h i h i .∗ (NumI : ( −1) :0 ) + mid h i .∗ ( 0 : ( + 1 ) : NumI ) ) . / NumI ;HiVI = f l i p l r ( HiVI ) ;% Now f i r s t the lower par t , we v e r i f y t ha t i n each v i n t e r v a l ,106% c w i l l shr ink , v w i l l increasef o r k = 1: leng th ( LoVI)−1v0 = [ LoVI ( k ) LoVI ( k + 1 ) ] ;z = Run pos (p , ds , v0 , theta ,0);%−ds posz1 = Run neg (p , ds , v0,− f l i p l r ( t he ta ) ,0) ;%+ds negi f ( z ==0 | | z1==0)data . pass = 0;d isp ( ’ f a i l i n lower par t , i n t e r v a l v : ’ )d isp ( v0 )endend% fo r the h igher par t , we v e r i f y t ha t i n each v i n t e r v a l ,% c w i l l shr ink , v w i l l decreasef o r k = 1: leng th ( HiVI )−1v0 = [ HiVI ( k ) HiVI ( k + 1 ) ] ;z = Run pos (p , ds , v0 , theta ,1);%−ds posz1 = Run neg (p , ds , v0,− f l i p l r ( t he ta ) ,1) ;%+ds negi f ( z ==0 | | z1==0)data . pass = 0;d isp ( ’ f a i l i n h igher par t , i n t e r v a l v : ’ )d isp ( v0 )endendi f ( data . pass==1)d isp ( ’ succeed in step four ’ ) ;e lsed isp ( ’ f a i l i n step four ’ ) ;endend% t h i s f unc t i on run i n mode1 u n t i l a t r a n s i t i o n to mode 2 ( the ta crosses 0% f a l l i n g ) or a t r a n s i t i o n to mode 4 ( c reaches c max ) or mode 3% ( c reaches c min ) t h i s more prec ise .% the ” or ” statement . s a t i s f a c t i o n o f one of the stop cond i t i on does not% garantee t ha t the i n i t a l ph i s empty .%% i t handle the t r a n s i t i o n manually to get the next s t a r t i n g pro jec tgon% inpu t : i n i t i a l pro jectgon ,% output : a c e l l data1 , con ta in ing the r e s u l t p ro jec tgon f o r mode2 ,% mode3 and mode4func t i on data1 = Run mode1 ( i n i tPh , invs , p )dim = 3; planes = [ 1 , 2 ; 1 , 3 ; 2 , 3 ] ;107i f ( nargin<1) % set de f au l t s f o r PLL model parametersopt . i n i t b b o x = [0 .90 0 .92 ;2 .3 2 .5 ;0 1 ] ;i n i t P h = ph createByBox ( dim , planes , opt . i n i t b b o x ) ;i n i t P h = ph conver t ( i n i tPh , ’ convex ’ ) ;endi f ( nargin<2)invs = [ ] ;endi f ( nargin<3)p = getP ( ) ;e lsep = getP ( p ) ;end%d isp lay the i n i t P h bounding boxdisp ( ’ en ter mode1 wi th [ c , v , the ta ] ’ ) ;d isp ( i n i t P h . bbox ) ;% resu l t , mode1 can t r a n s f e r to mode2( theta<0) and% mode4( theta>0 && c>cmax ) and mode3( theta<0 && c<cmin )data1 = c e l l ( 1 , 0 ) ;% set stop cond i t i onsca l lBacks . s l iceCond = @( i n f o ) ( f a l s e ) ; % we handle mode t r a n s i t i o n s manuallyca l lBacks . exi tCond = @( i n f o ) ex mode1 exitCond ( in fo , p ) ;% the c a l l to mode1 below computes the reachable sets u n t i l t he ta < 0ha = mode1 ( [ ] , ca l lBacks , i n i tPh , invs , p ) ;%load ( da ta f i l e name ) ;da ta f i l e name = [ ha . name, ’ ’ , ha . snames{1} , ’ . mat ’ ] ;[ reachData , t ime ] = LoadData ( da ta f i l e name ) ;i f ( t ime >= p .maxCompT)data1 . stop = 0;data1 . phin i t mode2 = [ ] ;data1 . phin i t mode3 = [ ] ;data1 . phin i t mode4 = [ ] ;d isp ( ’ compuation t ime exceeds the maxCompT, i t may not stop ’ )e lsedata1 . stop = 1;phsect1= ph createByBox ( dim , planes , . . .[ p . cmin , p . cmax ; p . vmin , p . vmax ; 0.0 ,2∗ p i ] ) ; % the ta >= 0phsect next2=ph createByBox ( dim , planes , . . .[ p . cmin , p . cmax ; p . vmin , p . vmax;−2∗pi , 0 . 0 ] ) ; % the ta <= 0phsect next3=ph createByBox ( dim , planes , . . .[ 0 , p . cmin ; p . vmin , p . vmax;−2∗pi , 0 ] ) ;% sat low108phsect next4=ph createByBox ( dim , planes , . . .[ p . cmax , p . cmax+0.1 ; p . vmin , p . vmax;0 ,2∗ p i ] ) ; % sat high%co l l e c t the r e s u l t pro jectgon , you may want i t as next s t a r t i n g%pro jec tagonph i n i t mode2 f r e f = [ ] ;ph i n i t mode3 f r e f = [ ] ;ph i n i t mode4 f r e f = [ ] ;% d i v i de the reachable sets i n t o two par ts :% the pa r t w i th the ta >= 0 , s to re i n curdata and data1{k , :}% the pa r t w i th the ta < 0 , compute the union o f these reg ions and use% tha t as the s t a r t i n g set f o r the next pa r t o f the ana lys i s .j = 1 ;f o r i = 1 : leng th ( reachData )i f ( ˜ isempty ( reachData{ i } ) )temp = ph canon ( ph i n t e r sec t ({ reachData{ i } ; phsect1 } ) ) ;i f ( ˜ isempty ( temp ) )data1 .mode1{ j }=temp ;j = j +1;endph cut 2 = ph canon ( ph i n t e r sec t ({ phsect next2 ; reachData{ i }} ) ) ;ph cut 3 = ph canon ( ph i n t e r sec t ({ phsect next3 ; reachData{ i }} ) ) ;ph cut 4 = ph canon ( ph i n t e r sec t ({ phsect next4 ; reachData{ i }} ) ) ;ph i n i t mode2 f r e f =ph union ({ ph i n i t mode2 f r e f ; ph cut 2 } ) ;%sat lowph i n i t mode3 f r e f =ph union ({ ph i n i t mode3 f r e f ; ph cut 3 } ) ;%sat highph i n i t mode4 f r e f =ph union ({ ph i n i t mode4 f r e f ; ph cut 4 } ) ;endend% p l o t% phs d isp lay ( data1 .mode1 , [ ] , [ ] , [ ] , ’ r ’ ) ; drawnow ;% get datadata1 . phin i t mode2 = ph i n i t mode2 f r e f ;data1 . phin i t mode3 = ph i n i t mode3 f r e f ;data1 . phin i t mode4 = ph i n i t mode4 f r e f ;endendfunc t i on exitCond = ex mode1 exitCond ( in fo , p )%de le te ph conta in ( i n f o . prevPh , i n f o . ph )exitCond= ph isempty ( i n f o . ph ) | | i n f o . ph . bbox ( 3 , 2 ) <0 | | . . .a l l ( d i f f ( i n f o . ph . bbox , [ ] ,2 )<1e−6) | | i n f o . compT>=p .maxCompT;109end % ex model exi tCond% t h i s f unc t i on run i n mode 2 u n t i l a t r a n s i t i o n to mode 1 ( the ta crosses 0% f a l l i n g ) or a t r a n s i t i o n to mode 4 ( c reaches c max ) or mode 3 ( c reaches% c min ) t h i s more prec ise .%% i t handle the t r a n s i t i o n manually to get the next s t a r t i n g pro jec tgon% inpu t : i n i t i a l pro jectgon ,% output : a c e l l data2 , con ta in ing the r e s u l t p ro jec tgon f o r% mode1 , mode3 and mode4func t i on data2 = Run mode2 ( i n i tPh , invs , p )dim = 3; planes = [ 1 , 2 ; 1 , 3 ; 2 , 3 ] ;i f ( nargin<1) % set de f au l t s f o r PLL model parametersopt . i n i t b b o x = [0.9075 0.9399;1.7803 1.7913;−0.0980 0 ] ;i n i t P h = ph createByBox ( dim , planes , opt . i n i t b b o x ) ;i n i t P h = ph conver t ( i n i tPh , ’ convex ’ ) ;endi f ( nargin<2)invs = [ ] ;endi f ( nargin<3)p = getP ( ) ;e lsep = getP ( p ) ;end%d isp lay the i n i t P h bounding boxdisp ( ’ en ter mode2 wi th [ c , v , the ta ] ’ ) ;d isp ( i n i t P h . bbox ) ;data2 = c e l l ( 1 , 0 ) ; % resu l t , mode2 can t r a n s f e r to mode1( theta >0) ,%mode3( theta<0 && c<c min ) and mode4( theta>0 && c>c max )% set e x i t cond i t i onsca l lBacks . s l iceCond = @( i n f o ) ( f a l s e );% we handle mode t r a n s i t i o n s manuallyca l lBacks . exi tCond = @( i n f o ) ex mode2 exitCond ( in fo , p ) ;% the c a l l to mode1 below computes the reachable sets u n t i l t he ta < 0ha = mode2 ( [ ] , ca l lBacks , i n i tPh , invs , p ) ; % the l a s t parameter means we% load datada ta f i l e name = [ ha . name, ’ ’ , ha . snames{1} , ’ . mat ’ ] ;[ reachData , t ime ] = LoadData ( da ta f i l e name ) ;i f ( t ime >= p .maxCompT)110data2 . stop = 0;data2 . phin i t mode1 = [ ] ;data2 . phin i t mode3 = [ ] ;data2 . phin i t mode4 = [ ] ;d isp ( ’ compuation t ime exceeds the maxCompT, i t may not stop ’ )e lsedata2 . stop = 1;phsect1= ph createByBox ( dim , planes , . . .[ 0 . 9 , 1 . 1 ; p . vmin , p . vmax;−2∗pi , 0 . 0 ] ) ; % the ta <= 0phsect next1=ph createByBox ( dim , planes , . . .[ 0 . 9 , 1 . 1 ; p . vmin , p . vmax ;0 .0 ,2∗ p i ] ) ; % the ta >= 0phsect next3=ph createByBox ( dim , planes , . . .[ 0 , 0 . 9 ; p . vmin , p . vmax;−2∗pi , 0 ] ) ; % the ta <= 0 & c<0.9phsect next4=ph createByBox ( dim , planes , . . .[ 1 . 1 , 1 . 2 ; p . vmin , p . vmax;0 ,2∗ p i ] ) ; % sat high% t h i s i s used to c o l l e c t the pro jec tgon f o r next 1 / f r e f cyc leph i n i t mode1 f r e f = [ ] ;ph i n i t mode3 f r e f = [ ] ;ph i n i t mode4 f r e f = [ ] ;% d i v i de the reachable sets i n t o two par ts :% the pa r t w i th the ta >= 0 , s to re i n curdata and data1{k , :}% the pa r t w i th the ta < 0 , compute the union o f these reg ions and use% tha t as the s t a r t i n g set f o r the next pa r t o f the ana lys i s .j = 1 ;f o r i = 1 : leng th ( reachData )i f ( ˜ isempty ( reachData{ i } ) )temp = ph canon ( ph i n t e r sec t ({ reachData{ i } ; phsect1 } ) ) ;i f ( ˜ isempty ( temp ) )data2 .mode2{ j }=temp ;j = j +1;endph cut 1 = ph canon ( ph i n t e r sec t ({ phsect next1 ; reachData{ i }} ) ) ;ph cut 3 = ph canon ( ph i n t e r sec t ({ phsect next3 ; reachData{ i }} ) ) ;ph cut 4 = ph canon ( ph i n t e r sec t ({ phsect next4 ; reachData{ i }} ) ) ;ph i n i t mode1 f r e f =ph union ({ ph i n i t mode1 f r e f ; ph cut 1 } ) ;%sat low modeph i n i t mode3 f r e f =ph union ({ ph i n i t mode3 f r e f ; ph cut 3 } ) ;%sat high modeph i n i t mode4 f r e f =ph union ({ ph i n i t mode4 f r e f ; ph cut 4 } ) ;endend% p l o t% phs d isp lay ( data2 .mode2 , [ ] , [ ] , [ ] , ’ b ’ ) ; drawnow ;% get some datadata2 . phin i t mode1 = ph i n i t mode1 f r e f ;data2 . phin i t mode3 = ph i n i t mode3 f r e f ;data2 . phin i t mode4 = ph i n i t mode4 f r e f ;111endendfunc t i on exitCond = ex mode2 exitCond ( in fo , p )exi tCond =( ph isempty ( i n f o . ph ) | | i n f o . ph . bbox ( 3 , 1 ) >0 | | . . .a l l ( d i f f ( i n f o . ph . bbox , [ ] ,2 )<1e−6)) | | i n f o . compT>=p .maxCompT;end% t h i s f unc t i on run mode3% mode3 i s the sat low mode , i t i s l i n e a r% inpu t : i n i t i a l pro jectgon ,% output : a c e l l data3 , con ta in ing the r e s u l t p ro jec tgon% f o r mode2 and mode4func t i on data3 = Run mode3 (p , i n i t P h )dim = 2; planes = [ 1 , 2 ] ;i f ( nargin<1)p = getP ( ) ;e lsep = getP ( p ) ;endi f ( nargin<2) % set de f au l t s f o r PLL model parametersopt . i n i t b b o x = [ p . vmin p . cmin∗p . f r e f ;−2∗p i 0 . 0 ] ;i n i t P h = ph createByBox ( dim , planes , opt . i n i t b b o x ) ;i n i t P h = ph conver t ( i n i tPh , ’ convex ’ ) ;end%d isp lay the i n i t P h bounding boxdisp ( ’ en ter mode1 wi th [ v , the ta ] ’ ) ;d isp ( i n i t P h . bbox ) ;data3 = c e l l ( 1 , 0 ) ;% set e x i t cond i t i onsca l lBacks . s l iceCond = @( i n f o ) ( f a l s e ) ; % we handle mode t r a n s i t i o n s manuallyca l lBacks . exi tCond = @( i n f o ) ex mode3 exitCond ( in fo , p ) ;ha = mode3 ( [ ] , ca l lBacks , i n i tPh , p , [ ] ) ;%c o l l e c t data from mode1da ta f i l e name = [ ha . name, ’ ’ , ha . snames{1} , ’ . mat ’ ] ;[ reachData , t ime ] = LoadData ( da ta f i l e name ) ;i f ( t ime >= p .maxCompT)112data3 . leavewa l l = 0 ;data3 . phin i t mode1 = [ ] ;d isp ( ’ compuation t ime exceeds the maxCompT, i t may not stop ’ )e lsedata3 . leavewa l l = 1 ;%the exitCond t e l l us t ha t i t w i l l leave the wa l l once stopphsect1= ph createByBox ( dim , planes , . . .[ p . vmin , p . vmax;−2∗pi , 0 . 0 ] ) ; % the ta <= 0phsect next1=ph createByBox ( dim , planes , . . .[ p . vmin , p . vmax ;0 .0 ,2∗ p i ] ) ; % the ta >= 0% t h i s i s used to c o l l e c t the pro jec tagon f o r next modeph i n i t mode1 f r e f = [ ] ;% d i v i de the reachable sets i n t o two par ts :% the pa r t w i th the ta >= 0 , s to re data3{k , :}% the pa r t w i th the ta < 0 , compute the union o f these reg ions and use% tha t as the s t a r t i n g set f o r the next pa r t o f the ana lys i s .j = 1 ;f o r i = 1 : leng th ( reachData )i f ( ˜ isempty ( reachData{ i } ) )temp = ph canon ( ph i n t e r sec t ({ reachData{ i } ; phsect1 } ) ) ;i f ( ˜ isempty ( temp ) )data3 .mode3{ j }=temp ;j = j +1;endph cut 1 = ph canon ( ph i n t e r sec t ({ phsect next1 ; reachData{ i }} ) ) ;ph i n i t mode1 f r e f =ph union ({ ph i n i t mode1 f r e f ; ph cut 1 } ) ;endend% p lo t , i t w i l l p l o t i n the x1 / x2 space , so I commented out% phs d isp lay ( data3 .mode3 , [ ] , [ ] , [ ] , ’ r ’ ) ; drawnow ;% get some datai f ( ˜ isempty ( ph i n i t mode1 f r e f ) )dim = 3; planes = [ 1 , 2 ; 1 , 3 ; 2 , 3 ] ;mode1 . i n i t b b o x = [ p . cmin p . cmin+1e−6 ; . . .ph i n i t mode1 f r e f . bbox ];% recover i t backin i t mode1 = ph createByBox ( dim , planes ,mode1 . i n i t b b o x ) ;in i t mode1 = ph conver t ( in i t mode1 , ’ convex ’ ) ;data3 . phin i t mode1 = in i t mode1 ;e lsedata3 . phin i t mode1 = [ ] ;enddisp ( ’ r e s u l t f o r next mode1 ’ )113disp ( data3 . phin i t mode1 ) ;endendfunc t i on exitCond = ex mode3 exitCond ( in fo , p )exi tCond = ph isempty ( i n f o . ph ) | | i n f o . ph . bbox ( 2 , 1 ) >0 | | . . .a l l ( d i f f ( i n f o . ph . bbox , [ ] ,2 )<1e−6) | | i n f o . compT>=p .maxCompT;end% t h i s f unc t i on run mode% mode4 i s the sat high mode , i t i s l i n e a r% inpu t : i n i t i a l pro jectagon ,% output : a c e l l data4 , con ta in ing the r e s u l t pro jec tagon f o r% mode2 and mode4func t i on data4 = Run mode4 (p , i n i t P h )dim = 2; planes = [ 1 , 2 ] ;i f ( nargin<1)p = getP ( ) ;e lsep = getP ( p ) ;endi f ( nargin<2) % set de f au l t s f o r PLL model parametersopt . i n i t b b o x = [ p . cmax∗p . f r e f p . vmax ;0 2∗p i ] ;i n i t P h = ph createByBox ( dim , planes , opt . i n i t b b o x ) ;i n i t P h = ph conver t ( i n i tPh , ’ convex ’ ) ;end%d isp lay the i n i t P h bounding boxdisp ( ’ en ter mode1 wi th [ c , v , the ta ] ’ ) ;d isp ( i n i t P h . bbox ) ;data4 = c e l l ( 1 , 0 ) ;% set e x i t cond i t i onsca l lBacks . s l iceCond = @( i n f o ) ( f a l s e ) ; % we handle mode t r a n s i t i o n s manuallyca l lBacks . exi tCond = @( i n f o ) ex mode4 exitCond ( in fo , p ) ;% the c a l l to mode1 below computes the reachable sets u n t i l t he ta < 0ha = mode4 ( [ ] , ca l lBacks , i n i tPh , p , [ ] ) ; % the l a s t parameter means we114%co l l e c t data from mode1da ta f i l e name = [ ha . name, ’ ’ , ha . snames{1} , ’ . mat ’ ] ;[ reachData , t ime ] = LoadData ( da ta f i l e name ) ;i f ( t ime >= p .maxCompT)data4 . leavewa l l = 0 ;data4 . phin i t mode2 = [ ] ;d isp ( ’ compuation t ime exceeds the maxCompT, i t may not stop ’ )e lse%the exitCond t e l l us t ha t i t w i l l leave the wa l l once stopdata4 . leavewa l l = 1 ;% the ta >= 0phsect1= ph createByBox ( dim , planes , [ p . vmin , p . vmax ;0 .0 ,2∗ p i ] ) ;% the ta <= 0phsect next2=ph createByBox ( dim , planes , [ p . vmin , p . vmax;−2∗pi , 0 . 0 ] ) ;% t h i s i s used to c o l l e c t the pro jec tgon f o r the next modeph i n i t mode2 f r e f = [ ] ;% d i v i de the reachable sets i n t o two par ts :% the pa r t w i th the ta >= 0 , s to re i n curdata and data4{k , :}% the pa r t w i th the ta < 0 , compute the union o f these reg ions and use% tha t as the s t a r t i n g set f o r the next pa r t o f the ana lys i s .j = 1 ;f o r i = 1 : leng th ( reachData )i f ( ˜ isempty ( reachData{ i } ) )temp = ph canon ( ph i n t e r sec t ({ reachData{ i } ; phsect1 } ) ) ;i f ( ˜ isempty ( temp ) )data4 .mode4{ j }=temp ; % i s data4 .mode3{ j } unneeded as we l l ?j = j +1;endph cut 2 = ph canon ( ph i n t e r sec t ({ phsect next2 ; reachData{ i }} ) ) ;ph i n i t mode2 f r e f =ph union ({ ph i n i t mode2 f r e f ; ph cut 2 } ) ;endend% p l o t% phs d isp lay ( data4 .mode4 , [ ] , [ ] , [ ] , ’ r ’ ) ; drawnow ;% get some datai f ( ˜ isempty ( ph i n i t mode2 f r e f ) )dim = 3; planes = [ 1 , 2 ; 1 , 3 ; 2 , 3 ] ;%recover i tmode2 . i n i t b b o x = [ p . cmax − 1e−6 p . cmax ; ph i n i t mode2 f r e f . bbox ] ;in i t mode2 = ph createByBox ( dim , planes ,mode2 . i n i t b b o x ) ;in i t mode2 = ph conver t ( in i t mode2 , ’ convex ’ ) ;data4 . phin i t mode2 = in i t mode2 ;115elsedata4 . phin i t mode2 = [ ] ;endendendfunc t i on exitCond = ex mode4 exitCond ( in fo , p )exi tCond = ph isempty ( i n f o . ph ) | | i n f o . ph . bbox (2 ,2)<0 | | . . .a l l ( d i f f ( i n f o . ph . bbox , [ ] ,2 )<1e−6) | | i n f o . compT>=p .maxCompT;end% posfunc t i on z=Run pos (p , dc , v0 , theta , p t )i f ( nargin<1)p = getP ( ) ;e lsep = getP ( p ) ;endi f ( nargin<2)dc = 0 .03 ;endi f ( nargin<3)%v0 = [1 . 8 1 . 8 2 ] ;%v0 = p . f r e f ∗ [ p . cmin−dc , p . cmin+dc ] ;%v0 = p . f r e f ∗ [ p . cmin−dc , p . cmax+dc ];%v0 = [2 .14 2 . 2 ] ;v0 = [2 .22 , 2 . 2 6 ] ;endi f ( nargin<4)the ta = [0 0 . 5 ] ;endi f ( nargin<5)%pt = 0;%lower pa r tp t = 1;%high pa r tend%%%% set the i n i t i a l reg ion%%%%%%%%opt . i n i t b b o x = [−dc ,0 ;0 ,1e−6; the ta ] ;% v i s very t inny , so t ha t i t i s c lose to zerosca l lBacks . exi tCond = @( i n f o ) ex mode1 exitCond ( i n f o ) ;ca l lBacks . s l iceCond = @( i n f o ) ( 0 ) ;dim = 3; planes = [ 1 , 2 ; 1 , 3 ; 2 , 3 ] ;i n i t P h = ph createByBox ( dim , planes , opt . i n i t b b o x ) ;i n i t P h = ph conver t ( i n i tPh , ’ convex ’ ) ;116%inv = e x p l l i n v ; invs = inv {1}; % the ta belongs to [0 , 2 p i ]invs = [ ] ;%%%% run the model and load data%%%%%ha = ex p l l 2 ( opt , v0 , cal lBacks , i n i tPh , invs , 1 ) ; %pos modeda ta f i l e name = [ ha . name, ’ ’ , ha . snames{1} , ’ . mat ’ ] ;reachData = LoadData ( da ta f i l e name ) ;ph sec t s tep1 = ph createByBox ( dim , planes , . . .[−1 ,1;−1 ,1;0 ,2∗ p i ] ) ;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%now get the next s t a r t i n g pro jec tagon :i n i t P h = [ ] ;r e s u l t s = c e l l (1 ,1) ;% c o l l e c t f o r p l o t s ;phsect=ph createByBox ( dim , planes ,[−1 ,1;−1 ,1;−2∗ pi , 0 ] ) ;%i n t e r s e c t w i th t h i s ph to get the next s t a r t i n g phf o r i = 1 : leng th ( reachData )%c o l l e c t r e s u l t s from har e su l t s{ i }=ph i n t e r sec t ({ reachData{ i } ; ph sec t s tep1 } ) ;%get the next s t a r t i n g phphu = ph i n t e r sec t ({ phsect ; reachData{ i }} ) ;i n i t P h = ph union ({phu ; i n i t P h } ) ;end% containment proper ty , check i f the f i n a l% r e s u l t conta ined i n the checkboxcheckbox = ph createByBox ( dim , planes , . . .[−dc , dc;−1,1;− f l i p l r ( t he ta ) ] ) ;y1 = ph conta in ( checkbox , i n i t P h ) ;i f ( p t ==0)y2 = i n i t P h . bbox(2 ,1)>0;% lower par t , v increaseelsey2 = i n i t P h . bbox (2 ,2)<0; % higher par t , v decrease ,endi f ( y1 && y2)%the second cond i t i on shows t ha t the f low w i l l go updisp ( ’ s a t i s f i e d ’ ) ;z=1;e lsez=0;d isp ( ’ no ’ ) ;endendfunc t i on exitCond = ex mode1 exitCond ( i n f o )exi tCond = ph isempty ( i n f o . ph ) | | i n f o . ph . bbox ( 3 , 1 ) <0 | | . . .117a l l ( d i f f ( i n f o . ph . bbox , [ ] ,2 )<1e−6);endfunc t i on invs = e x p l l i n vA1 = [0 ,0 ,−1;0 ,0 ,1 ] ;A2 = [0 ,0 ,1 ;0 ,0 ,−1] ;b1 = [0 ;2∗ p i ] ;b2 = [0 ;2∗ p i ] ;inv1 = l p c r ea t e (A1 , b1);%posinv2 = l p c r ea t e (A2 , b2);%neginvs = { inv1 ; inv2 } ;end%negfunc t i on z=Run neg (p , dc , v0 , theta , p t )i f ( nargin<1)p = getP ( ) ;e lsep = getP ( p ) ;endi f ( nargin<2)dc = 0 .03 ;endi f ( nargin<3)%v0 = [1 . 8 1 . 8 2 ] ;%v0 = p . f r e f ∗ [ p . cmin−dc , p . cmin+dc ] ;%v0 = p . f r e f ∗ [ p . cmax−dc , p . cmax+dc ] ;v0 =[2 .22 2 . 2 6 ] ;endi f ( nargin<4)the ta = [−0.5 0 ] ;endi f ( nargin<5)%pt = 0;%lower pa r tp t = 1;%high pa r tend% v i s very t inny , so t ha t i t i s c lose to zerosopt . i n i t b b o x = [0 , dc ;0 ,1e−6; the ta ] ;ca l lBacks . exi tCond = @( i n f o ) ex mode2 exitCond ( i n f o ) ;ca l lBacks . s l iceCond = @( i n f o ) ( 0 ) ;dim = 3; planes = [ 1 , 2 ; 1 , 3 ; 2 , 3 ] ;i n i t P h = ph createByBox ( dim , planes , opt . i n i t b b o x ) ;i n i t P h = ph conver t ( i n i tPh , ’ convex ’ ) ;%inv = e x p l l i n v ; invs = inv {2};i nvs = [ ] ;118ha = ex p l l 2 ( opt , v0 , cal lBacks , i n i tPh , invs , 2 ) ; % mode2da ta f i l e name = [ ha . name, ’ ’ , ha . snames{1} , ’ . mat ’ ] ;reachData = LoadData ( da ta f i l e name ) ;ph sec t s tep1 = ph createByBox ( dim , planes , . . .[−1,1;−1,1;−2∗pi , 0 ] ) ;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%now get the next s t a r t i n g pro jec tagon :i n i t P h = [ ] ;r e s u l t s = c e l l (1 ,1) ;% c o l l e c t f o r p l o t s ;phsect=ph createByBox ( dim , planes ,[−1 ,1;−1 ,1;0 ,2∗ p i ] ) ;%i n t e r s e c t w i th t h i s ph to get the next s t a r t i n g phf o r i = 1 : leng th ( reachData )%c o l l e c t r e s u l t s from har e su l t s{ i }=ph i n t e r sec t ({ reachData{ i } ; ph sec t s tep1 } ) ;%get the next s t a r t i n g phphu = ph i n t e r sec t ({ phsect ; reachData{ i }} ) ;i n i t P h = ph union ({phu ; i n i t P h } ) ;endcheckbox = ph createByBox ( dim , planes , . . .[−dc , dc;−1,1;− f l i p l r ( t he ta ) ] ) ;y1=ph conta in ( checkbox , i n i t P h ) ;i f ( p t ==0)y2 = i n i t P h . bbox(2 ,1)>0;% lower par t , v increaseelsey2 = i n i t P h . bbox (2 ,2)<0; % higher par t , v decrease ,endi f ( y1 && y2)%the second cond i t i on shows the f low w i l l go upz = 1;d isp ( ’ s a t i s f i e d ’ ) ;e lsez=0;d isp ( ’ no ’ ) ;end%cra c lose ;end%func t i on exitCond = ex mode2 exitCond ( i n f o )exi tCond = ph isempty ( i n f o . ph ) | | i n f o . ph . bbox ( 3 , 2 ) >0 | | . . .a l l ( d i f f ( i n f o . ph . bbox , [ ] ,2 )<1e−6);119endfunc t i on invs = e x p l l i n vA1 = [0 ,0 ,−1;0 ,0 ,1 ] ;A2 = [0 ,0 ,1 ;0 ,0 ,−1] ;b1 = [0 ;2∗ p i ] ;b2 = [0 ;2∗ p i ] ;inv1 = l p c r ea t e (A1 , b1 ) ; inv2 = l p c r ea t e (A2 , b2 ) ;invs = { inv1 ; inv2 } ;end%mode 1 , theta>0func t i on ha1 = mode1( opt , ca l lBacks , i n i tPh , invs , p )i f ( nargin<1) % use de f au l t valuesopt = [ ] ;ca l lBacks = [ ] ;i n i t P h = [ ] ;i nvs = [ ] ;p = getP ( p ) ;end % use de f au l t valueha1 = ex p l l ha1 ( opt , ca l lBacks , i n i tPh , invs , p ) ;ha1 = ha reach ( ha1 ) ; % run Coho u n t i l we swi tch to mode2end % mode1func t i on ha1 = ex p l l ha1 ( opt , ca l lBacks , i n i tPh , invs , p )i f ( nargin<1) opt = [ ] ; endi f ( isempty ( opt ) | | ˜ i s f i e l d ( opt , ’ i n i t bbox ’ ) | | . . .isempty ( opt . i n i t b b o x ) )opt . i n i t b b o x = [ 0 . 9 , 0 . 900001 ; 2 . 7 , 2 . 8 ; 0 , 5 ] ;endi f ( isempty ( opt ) | | ˜ i s f i e l d ( opt , ’ i n i t s t a t e ’ ) | | . . .isempty ( opt . i n i t s t a t e ) )opt . i n i t s t a t e = 1;endi f ( narg in < 2) ca l lBacks = [ ] ; end ;i f ( ( narg in < 3) | | isempty ( i n i t P h ) )e r r o r ( ’ i n i t P h must be provided ’ ) ;endi f ( nargin<4) invs = [ ] ; endi f ( nargin<5)p = getP ( p ) ;endphOpt = e x p l l o p t ;s ta tes ( 1 ) = ha s ta te ( ’mode1 ’ ,@( l p ) ( ex p l l mode l ( p , l p ) ) , . . .invs , phOpt , ca l lBacks ) ;120t r ans = [ ] ;source = ’mode1 ’ ;ha1 = ha create ( ’ p l l 1 ’ , s ta tes , t rans , source , i n i t P h ) ;end % ex p l l ha1f unc t i on l d i = ex p l l mode l ( p , l p )x=1; y=2; z=3;bbox = lp box ( l p ) ;avgs = mean( bbox , 2 ) ; ex t ras = d i f f ( bbox , [ ] , 2 ) / 2 ;% mode1 x(3)>=0k = −1;i f ( nargin<1)p = getP ( p ) ;enda = k∗p . g1 ; b = p . g2∗p . f r e f ; c = −b∗p . cc ; d = 1/p .N;e = −p .K ; f = −p . f r e f ;A1 = [0 , 0 , 0 ; . . .b , 0 , 0 ; . . .0 ,0 ,e ] ;B1 = [ a ; c ; f ] ;U1 = zeros ( 3 , 1 ) ;A2 = [0 , 0 , 0 ; . . .0 ,0 ,0 ; . . .−avgs ( y ) / avgs ( x ) ˆ 2 , 1 / avgs ( x ) , 0 ] ;B2 = [ 0 ; 0 ; avgs ( y ) / avgs ( x ) ] ;% e r r = ( avgs ( y ) dxˆ2−avgs ( x )∗dx∗dy ) / avgs ( x ) ˆ 2 ( avgs ( x )+dx )% dx=[−h1 , h1 ] ; dy=[−h2 ; h2 ] ;ep11 = ( avgs ( y )∗ ex t ras ( x)ˆ2−avgs ( x )∗ ex t ras ( x )∗ ex t ras ( y ) ) / . . .avgs ( x ) ˆ 2∗ ( avgs ( x)−ex t ras ( x ) ) ;ep12 = ( avgs ( y )∗ ex t ras ( x ) ˆ2+ avgs ( x )∗ ex t ras ( x )∗ ex t ras ( y ) ) / . . .avgs ( x ) ˆ 2∗ ( avgs ( x)−ex t ras ( x ) ) ;ep21 = ( avgs ( y )∗ ex t ras ( x ) ˆ2+ avgs ( x )∗ ex t ras ( x )∗ ex t ras ( y ) ) / . . .avgs ( x ) ˆ 2∗ ( avgs ( x )+ ex t ras ( x ) ) ;ep22 = ( avgs ( y )∗ ex t ras ( x)ˆ2−avgs ( x )∗ ex t ras ( x )∗ ex t ras ( y ) ) / . . .avgs ( x ) ˆ 2∗ ( avgs ( x )+ ex t ras ( x ) ) ;exx1 = −avgs ( x ) ˆ 2 / ( 4∗ avgs ( y ) )∗ ex t ras ( y ) ˆ 2 / . . .avgs ( x ) ˆ 2∗ ( avgs ( x)−ex t ras ( x ) ) ;exx2 = −avgs ( x ) ˆ 2 / ( 4∗ avgs ( y ) )∗ ex t ras ( y ) ˆ 2 / . . .avgs ( x ) ˆ 2∗ ( avgs ( x )+ ex t ras ( x ) ) ;eps = [ ep11 ; ep12 ; ep21 ; ep22 ; exx1 ; exx2 ] ;e r r = [ min ( eps ) ; max( eps ) ] ;B2( z ) = B2( z ) + mean( e r r ) ;e r r = e r r − mean( e r r ) ; % make e r r o r balanceU2 = [ 0 ; 0 ; abs ( e r r ( 1 ) ) ] ;121A = A1+d∗A2 ; B = B1+d∗B2 ; U = U1+d∗U2;U = U+1e−9; %only necessary when ob jec t i s face−none / he igh tl d i = i n t c r e a t e (A,B,U ) ;end % exp pl l mode1func t i on phOpt = e x p l l o p tphOpt . fwdOpt = ph getOpt ;phOpt . fwdOpt . ob jec t = ’ ph ’ ;% 1 model f o r whole ph , make the computat ion f a s t% phOpt . fwdOpt . model = ’ t imeStep ’ ; % slow down f o r d i sp lay% phOpt . fwdOpt . t imeStep = 0 . 1 ;end % ex p l l o p t%mode2 theta<0func t i on ha2 = mode2( opt , ca l lBacks , i n i tPh , inv , p )i f ( nargin<1) % use de f au l t valuesopt = [ ] ;ca l lBacks = [ ] ;i n i t P h = [ ] ;i nv = [ ] ;p = getP ( p ) ;end % use de f au l t valueha2 = ex p l l ha2 ( opt , ca l lBacks , i n i tPh , inv , p ) ;ha2 = ha reach ( ha2 ) ;endfunc t i on ha2 = ex p l l ha2 ( opt , ca l lBacks , i n i tPh , inv , p )i f ( nargin<1) , opt = [ ] ; endi f ( isempty ( opt ) | | ˜ i s f i e l d ( opt , ’ i n i t bbox ’ ) . . .| | isempty ( opt . i n i t b b o x ) )opt . i n i t b b o x = [ 0 . 9 , 0 . 9 001 ; 2 . 7 , 2 . 8 ; 0 , 5 ] ;endi f ( isempty ( opt ) | | ˜ i s f i e l d ( opt , ’ i n i t s t a t e ’ ) . . .| | isempty ( opt . i n i t s t a t e ) )opt . i n i t s t a t e = 1;endi f ( narg in < 2) ca l lBacks = [ ] ; end ;i f ( ( narg in < 3) | | isempty ( i n i t P h ) )e r r o r ( ’ i n i t P h must be provided ’ ) ;endi f ( nargin<4) inv = [ ] ; endi f ( nargin<5)p = getP ( p ) ;endphOpt = e x p l l o p t 2 ;122s ta tes ( 1 ) = ha s ta te ( ’mode2 ’ ,@( l p ) ( ex pl l mode2 (p , l p ) ) . . ., inv , phOpt , ca l lBacks ) ;t rans = [ ] ;source = ’mode2 ’ ;ha2 = ha create ( ’ p l l 1 ’ , s ta tes , t rans , source , i n i t P h ) ;endfunc t i on l d i = ex pl l mode2 (p , l p )x=1; y=2; z=3;bbox = lp box ( l p ) ;avgs = mean( bbox , 2 ) ; ex t ras = d i f f ( bbox , [ ] , 2 ) / 2 ;i f ( nargin<1)p = getP ( p ) ;end% mode2 x(3)<=0k = 1;a = k∗p . g1 ; b = p . g2∗p . f r e f ; c = −b∗p . cc ; d = 1/p .N;e = −p .K ; f = −p . f r e f ;A1 = [0 , 0 , 0 ; . . .b , 0 , 0 ; . . .0 ,0 ,e ] ;B1 = [ a ; c ; f ] ;U1 = zeros ( 3 , 1 ) ;A2 = [0 , 0 , 0 ; . . .0 ,0 ,0 ; . . .−avgs ( y ) / avgs ( x ) ˆ 2 , 1 / avgs ( x ) , 0 ] ;B2 = [ 0 ; 0 ; avgs ( y ) / avgs ( x ) ] ;% e r r = ( avgs ( y ) dxˆ2−avgs ( x )∗dx∗dy ) / avgs ( x ) ˆ 2 ( avgs ( x )+dx )% dx=[−h1 , h1 ] ; dy=[−h2 ; h2 ] ;ep11 = ( avgs ( y )∗ ex t ras ( x)ˆ2−avgs ( x )∗ ex t ras ( x )∗ ex t ras ( y ) ) / . . .avgs ( x ) ˆ 2∗ ( avgs ( x)−ex t ras ( x ) ) ;ep12 = ( avgs ( y )∗ ex t ras ( x ) ˆ2+ avgs ( x )∗ ex t ras ( x )∗ ex t ras ( y ) ) / . . .avgs ( x ) ˆ 2∗ ( avgs ( x)−ex t ras ( x ) ) ;ep21 = ( avgs ( y )∗ ex t ras ( x ) ˆ2+ avgs ( x )∗ ex t ras ( x )∗ ex t ras ( y ) ) / . . .avgs ( x ) ˆ 2∗ ( avgs ( x )+ ex t ras ( x ) ) ;ep22 = ( avgs ( y )∗ ex t ras ( x)ˆ2−avgs ( x )∗ ex t ras ( x )∗ ex t ras ( y ) ) / . . .avgs ( x ) ˆ 2∗ ( avgs ( x )+ ex t ras ( x ) ) ;exx1 = −avgs ( x ) ˆ 2 / ( 4∗ avgs ( y ) )∗ ex t ras ( y ) ˆ 2 / . . .avgs ( x ) ˆ 2∗ ( avgs ( x)−ex t ras ( x ) ) ;exx2 = −avgs ( x ) ˆ 2 / ( 4∗ avgs ( y ) )∗ ex t ras ( y ) ˆ 2 / . . .avgs ( x ) ˆ 2∗ ( avgs ( x )+ ex t ras ( x ) ) ;123eps = [ ep11 ; ep12 ; ep21 ; ep22 ; exx1 ; exx2 ] ;e r r = [ min ( eps ) ; max( eps ) ] ;B2( z ) = B2( z ) + mean( e r r ) ;e r r = e r r − mean( e r r );% make e r r o r balanceU2 = [ 0 ; 0 ; abs ( e r r ( 1 ) ) ] ;A = A1+d∗A2 ; B = B1+d∗B2 ; U = U1+d∗U2;U = U+1e−9; %only necessary when ob jec t i s face−none / he igh tl d i = i n t c r e a t e (A,B,U ) ;endfunc t i on phOpt = e x p l l o p t 2phOpt . fwdOpt = ph getOpt ;phOpt . fwdOpt . ob jec t = ’ ph ’ ;%make the computat ion f a s tend% t h i s i s the sa t u r a t i on mode , I f i r s t model the lower par t ,% which means when PLL% sat ruate , i t w i l l sa tu ra te a t the Cmin = 0 . 9 .f unc t i on ha1 = mode3( opt , ca l lBacks , i n i tPh , p , invs )i f ( nargin<1) % use de f au l t valuesopt = [ ] ;ca l lBacks = [ ] ;i n i t P h = [ ] ;p = getP ( p ) ;invs = [ ] ;end % use de f au l t valueha1 = ex p l l ha1 ( opt , ca l lBacks , i n i tPh , p , invs ) ;ha1 = ha reach ( ha1 ) ; % run Coho u n t i l we swi tch to mode1end % mode1func t i on ha1 = ex p l l ha1 ( opt , ca l lBacks , i n i tPh , p , invs )i f ( nargin<1) opt = [ ] ; endi f ( isempty ( opt ) | | ˜ i s f i e l d ( opt , ’ i n i t bbox ’ ) . . .| | isempty ( opt . i n i t b b o x ) )opt . i n i t b b o x = [ 2 . 7 , 2 . 8 ; 0 , 5 ] ;%opt . i n i t b b o x = [0 .9 ,1 .1 ;3 .0 ,3 .2 ; −2∗ pi , 0 ] ;endi f ( isempty ( opt ) | | ˜ i s f i e l d ( opt , ’ i n i t s t a t e ’ ) | | . . .isempty ( opt . i n i t s t a t e ) )opt . i n i t s t a t e = 1;endi f ( narg in < 2) ca l lBacks = [ ] ; end ;i f ( ( narg in < 3) | | isempty ( i n i t P h ) )e r r o r ( ’ i n i t P h must be provided ’ ) ;end124i f ( nargin<4)p = getP ( p ) ;endi f ( nargin<5)invs = [ ] ;endphOpt = e x p l l o p t ;s ta tes ( 1 ) = ha s ta te ( ’mode3 ’ ,@( l p ) ( ex pl l mode3 (p , l p ) ) , . . .invs , phOpt , ca l lBacks ) ;t rans = [ ] ;source = ’mode3 ’ ;ha1 = ha create ( ’ p l l 1 ’ , s ta tes , t rans , source , i n i t P h ) ;end % ex p l l ha1f unc t i on l d i = ex pl l mode3 (p , l p )% mode3 x(3)<=0k = 0;i f ( nargin<1)p = getP ( p ) ;enda = k ; b = p . g2∗p . f r e f ∗(p . cmin−p . cc ) ; s = 1 / p . cmin ;d = 1/p .N; e = −p .K ; f = −p . f r e f ;A = [ 0 , 0 ; . . .s , e ] ;B = [ b ; f ] ;%U = zeros ( 2 , 1 ) ;U = [1e−9;1e−9];l d i = i n t c r e a t e (A,B,U ) ;end % exp pl l mode1func t i on phOpt = e x p l l o p tphOpt . fwdOpt = ph getOpt ;phOpt . fwdOpt . ob jec t = ’ ph ’ ; %make the computat ion f a s t% phOpt . fwdOpt . model = ’ t imeStep ’ ; % slow down f o r d i sp lay% phOpt . fwdOpt . t imeStep = 0 . 1 ;end % ex p l l o p t% t h i s i s the sa t u r a t i on mode , model the h igher par t ,% which means when PLL% sat ruate , i t w i l l sa tu ra te a t the Cmax = 1 . .f unc t i on ha1 = mode4( opt , ca l lBacks , i n i tPh , p , invs )i f ( nargin<1) % use de f au l t values125opt = [ ] ;ca l lBacks = [ ] ;i n i t P h = [ ] ;p = getP ( p ) ;invs = [ ] ;end % use de f au l t valueha1 = ex p l l ha1 ( opt , ca l lBacks , i n i tPh , p , invs ) ;ha1 = ha reach ( ha1 ) ; % run Coho u n t i l we swi tch to mode1end % mode1func t i on ha1 = ex p l l ha1 ( opt , ca l lBacks , i n i tPh , p , invs )i f ( nargin<1) opt = [ ] ; endi f ( isempty ( opt ) | | ˜ i s f i e l d ( opt , ’ i n i t bbox ’ ) | | . . .isempty ( opt . i n i t b b o x ) )opt . i n i t b b o x = [ 2 . 7 , 2 . 8 ; 0 , 5 ] ;%opt . i n i t b b o x = [0 .9 ,1 .1 ;3 .0 ,3 .2 ; −2∗ pi , 0 ] ;endi f ( isempty ( opt ) | | ˜ i s f i e l d ( opt , ’ i n i t s t a t e ’ ) | | . . .isempty ( opt . i n i t s t a t e ) )opt . i n i t s t a t e = 1;endi f ( narg in < 2) ca l lBacks = [ ] ; end ;i f ( ( narg in < 3) | | isempty ( i n i t P h ) )e r r o r ( ’ i n i t P h must be provided ’ ) ;endi f ( nargin<4)p = getP ( p ) ;endi f ( nargin<5)invs = [ ] ;endphOpt = e x p l l o p t ;s ta tes ( 1 ) = ha s ta te ( ’mode4 ’ ,@( l p ) ( ex pl l mode4 (p , l p ) ) , . . .invs , phOpt , ca l lBacks ) ;t rans = [ ] ;source = ’mode4 ’ ;ha1 = ha create ( ’ p l l 1 ’ , s ta tes , t rans , source , i n i t P h ) ;end % ex p l l ha1f unc t i on l d i = ex pl l mode4 (p , l p )% mode4 x(3)>=0k = 0;i f ( nargin<1)p = getP ( p ) ;end126a = k ; b = p . g2∗p . f r e f ∗(p . cmax−p . cc ) ; c = 1 / p . cmax ;d = 1/p .N; e = −p .K ; f = −p . f r e f ;A = [0 ,0 ; . . .c , e ] ;B = [ b ; f ] ;%U = zeros ( 2 , 1 ) ;U = [1e−9;1e−9];l d i = i n t c r e a t e (A,B,U ) ;end % exp pl l mode1func t i on phOpt = e x p l l o p tphOpt . fwdOpt = ph getOpt ;phOpt . fwdOpt . ob jec t = ’ ph ’ ;% phOpt . fwdOpt . model = ’ t imeStep ’ ; % slow down f o r d i sp lay% phOpt . fwdOpt . t imeStep = 0 . 1 ;end % ex p l l o p t% the aim of t h i s f i l e i s to s imu la te the p l l but w i th no t r a n s i t i o n s and% make i t centered at ( c0 , v0 ) .% c = c0 + c1 ;% v = v0 + v1 ;% v1 (0 ) = 0 ; c1 (0 ) = 0% c1dot = s igh ( ph ) , v1dot = g2∗ f r e f ( c0+c1−ccode ) = g2∗ f r e f ( c1 − ( ccode−c0 ) )% phdot = 1 /N∗( v0+v1 ) / ( c0+c1 ) − f r e f − k∗phfunc t i on ha = ex p l l 2 ( opt , v0 , cal lBacks , i n i tPh , invs ,mode)i f ( nargin<1) , opt = [ ] ; end % use de f au l t valueha = ex p l l h a ( opt , v0 , cal lBacks , i n i tPh , invs ,mode ) ;ha = ha reach ( ha ) ;endfunc t i on ha = ex p l l h a ( opt , v0 , cal lBacks , i n i tPh , invs ,mode)i f ( nargin<1) , opt = [ ] ; endi f ( isempty ( opt ) | | ˜ i s f i e l d ( opt , ’ i n i t bbox ’ ) | | isempty ( opt . i n i t b b o x ) )opt . i n i t b b o x = [ 0 , 0 ; 0 , 0 ; 0 , 0 . 0 1 ] ;endi f ( isempty ( opt ) | | ˜ i s f i e l d ( opt , ’ i n i t s t a t e ’ ) | | isempty ( opt . i n i t s t a t e ) )opt . i n i t s t a t e = 1;endi f ( nargin<2)v0 =[p . f r e f ∗(p . cmin−0.01) p . f r e f ∗(p . cmax+0 . 01 ) ] ;endi f ( nargin<3)ca l lBacks . s l iceCond = @( i n f o ) ( 0 ) ;end127i f ( nargin<4)e r r o r ( ’ i n i t i a l pro jec tagon cannot be empty ! ’ ) ;endi f ( nargin<5)invs = [ ] ;endi f ( nargin<6)mode = 1;endphOpt = e x p l l o p t ;i f (mode==1)s ta tes = ha s ta te ( ’mode1 ’ ,@( l p ) ( ex p l l mode l ( lp ,mode , v0 ) ) . . ., invs , phOpt , ca l lBacks ) ;source = ’mode1 ’ ;e lses ta tes = ha s ta te ( ’mode2 ’ ,@( l p ) ( ex p l l mode l ( lp ,mode , v0 ) ) , . . .invs , phOpt , ca l lBacks ) ;source = ’mode2 ’ ;endt rans = [ ] ;ha = ha create ( ’ p l l ’ , s ta tes , t rans , source , i n i t P h ) ;end%%%please look a t mark ’ s model ( s t r a t egy . pdf and example2 . pdf )f unc t i on l d i = ex p l l mode l ( lp ,mode , v0 )x=1; y=2; z=3;swi tch (mode)case 1 % x(3)>=0s = −1;case 2 % x(3)<=0s = 1;endp = getP ( ) ;c0= v0 / p . f r e f ;cmid = mean( c0 ) ; vmid = mean( v0 ) ; c t = d i f f ( c0 ) / ( 2∗ cmid ) ;v t = d i f f ( v0 ) / ( 2∗ vmid ) ; c t2 = 1−c t ˆ 2 ;e1 =c t∗ v t ; e2 = c t ˆ 2 / ( 2∗ c t2 ) ; a = 1+e2 ;e3 = e1+(1+ v t +c t +2∗e1)∗e2 ;A = [0 , 0 , 0 ; . . .p . g2∗p . f r e f , 0 , 0 ; . . .−(a∗p . f r e f ) / cmid , ( a∗p . f r e f ) / vmid , −p .K ] ;B = [ s∗p . g1 ; −p . f r e f ∗p . g2∗(p . cc − vmid / p . f r e f ) ; ( a−1)∗p . f r e f ] ;U = [ 0 ; abs ( p . g2∗ v t∗vmid ) ; e3∗p . f r e f ] ;l d i = i n t c r e a t e (A,B,U ) ;128endfunc t i on phOpt = e x p l l o p tphOpt . fwdOpt = ph getOpt ;% 1 model f o r whole ph , make the computat ion f a s tphOpt . fwdOpt . ob jec t = ’ ph ’ ;% phOpt . fwdOpt . model = ’ t imeStep ’ ; % slow down f o r d i sp lay% phOpt . fwdOpt . t imeStep = 0 .05 ;endfunc t i on varargout = get params (p , va ra rg in )i n = va ra rg in ;out = c e l l ( l eng th ( i n ) / 2 , 1 ) ;f o r ( i = 1 : 2 : leng th ( i n ) )i f ( ˜ isempty ( p ) && i s f i e l d ( p , i n{ i } ) )out {( i +1)/2} = g e t f i e l d ( p , i n{ i } ) ;e lse out {( i +1)/2} = in{ i +1};endendvarargout = out ;endf unc t i on y = getP ( p )i f ( nargin<1)p = [ ] ;endy . f r e f = get params (p , ’ f r e f ’ , 2 ) ;y . g1 = get params (p , ’ g1 ’ ,−0.01) ;y . g2 = get params (p , ’ g2 ’ ,−0.002) ;y . cmin = get params (p , ’ cmin ’ , 0 . 9 ) ;y . cmax = get params (p , ’ cmax ’ , 1 . 1 ) ;y .K = get params (p , ’ K ’ , 0 . 8 ) ;y . vmin = get params (p , ’ vmin ’ , 1 . 5 ) ;y . vmax = get params (p , ’ vmax ’ , 2 . 5 ) ;y .N = get params (p , ’N’ , 1 . 0 ) ;y . cc = get params (p , ’ cc ’ , 1 . 0 ) ;y .maxCompT = get params (p , ’maxCompT’ , 5 0 ) ;end% t h i s f unc t i on i s used to load the data f o r% a l l Run modei ( ) f unc t i onsf unc t i on [ data , t ime ] = LoadData (name)i f ( nargin<1)e r r o r ( ’ to load the data , you must prov ide the name ’ ) ;end129load (name ) ;data = reachData . sets ;t ime = reachData . compT ; % computat ion t imeend% t h i s f unc t i on w i l l so lve a l i n e a r program to check whether :% s t a r t i n g from a lower s t r i p i t w i l l move upward or% s t a r t i n g from a higher s t r i p i t w i l l move downward% inpu t : ph the progjecgon to check% f r e f l ock ing frequency% a=0 lower s t r i p , a=1 h igher s t r i p% output : ds , the minimum and maximum dis tance to the l i n e v / c= f r e ff unc t i on ds=checkPass ( ph , f r e f , a )% v = f r e f ∗( c−ds ) we are i n t e r es t ed i n the min and max value o f ds% so l e t i t be the t a r ge t f unc t i on% ds = c − v / f r e f , when a = 0% Likewise , when a = 1 , ds = −c +v / f r eq ;i f ( nargin<1)ph = [ ] ;endi f ( nargin<2)f r e f = 2 ;endi f ( nargin<3)a = 0;endA = ph . hu l lLP .A ;b = ph . hu l lLP . b ;Aeq = ph . hu l lLP . Aeq ;beq = ph . hu l lLP . beq ;i f ( a==0)f1 = [1 , −1/ f r e f , 0 ] ;f2 = [−1 ,1/ f r e f , 0 ] ;e lsef1 = [−1 ,1/ f r e f , 0 ] ;f2 = [1 ,−1/ f r e f , 0 ] ;end[ ˜ , ds1 ] = l i n p r og ( f1 ,A, b , Aeq , beq ) ;[ ˜ , ds2 ] = l i n p r og ( f2 ,A, b , Aeq , beq ) ;ds = [ ds1,−ds2 ] ;end130

Cite

Citation Scheme:

        

Citations by CSL (citeproc-js)

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}]}"
                            async >
                            </script>
                            </div>
                        
                    
IIIF logo Our image viewer uses the IIIF 2.0 standard. To load this item in other compatible viewers, use this url:
http://iiif.library.ubc.ca/presentation/dsp.24.1-0165974/manifest

Comment

Related Items