**Simulation of Piecewise Hybrid Dynamical Systems in Matlab**

Fatima El Guezar and Hassane Bouzahir

Additional information is available at the end of the chapter

http://dx.doi.org/10.5772/48570

## **1. Introduction**

A hybrid dynamical system is a system containing on the same time continuous state variables and event variables in interaction. We find hybrid systems in different fields. We cite robotic systems, chemical systems controlled by vans and pumps, biological systems (growth and division) and nonlinear electronics systems.

Because of interaction between continuous and discrete aspects, the behavior of hybrid systems can be seen as extremely complex. However, this behavior becomes relatively simple for piece-wise affine hybrid dynamical systems that can, in contrast, generate bifurcation and chaos. There are many examples such as power electronics DC-DC converters.

The common power electronics DC-DC converters are the buck converter and the boost converter. They are switching systems with time variant structure [9].

DC-DC converters are widely used in industrial, commercial, residential and aerospace environments. These circuits are typically controlled by PWM (Piece Wise Modulation) or other similar techniques to regulate the tension and the current given to the charges. The controller decides to pass from one configuration to another by considering that transitions occur cyclically or in discrete time. In order to make the analysis possible, most of mathematical treatments use some techniques that are based on averaging or discretization. Averaging can mean to wrong conclusions on operation of a system. Discrete models do not give any information on the state of the system between the sampled instants. In addition, they are difficult to obtain. In fact, in most cases, a pure analytic study is not possible. Another possible approach to analyze these converters can be done via some models of hybrid dynamical systems. DC-DC converters are particularly good candidates for this type of analysis because of their natural hybrid structure. The nature of commutations of these systems makes them strongly nonlinear. They present specific complex phenomena such as fractals structures of bifurcation and chaos.

©2012 El Guezar and Bouzahir, licensee InTech. This is an open access chapter distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited. ©2012 El Guezar and Bouzahir, licensee InTech. This is a paper distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

2 Will-be-set-by-IN-TECH 4 MATLAB – A Fundamental Tool for Scienti c Computing and Engineering Applications – Volume 3 Simulation of Piecewise Hybrid Dynamical Systems in Matlab <sup>3</sup>

The study of nonlinear dynamics of DC-DC converters started in 1984 by works of Brockett and Wood [4]. Since then, chaos and nonlinear dynamics in power electronics circuits have attracted different research groups around the world. Different nonlinear phenomena have been observed such as routes to chaos following the period doubling cascade [16], [5], [19], [20] and [23] or quasi-periodic phenomena [6], [7] and [8], besides border collision bifurcations [23] and [2].

present the algorithm that detects events' occurrence and devote a subsection to our approach efficiency. Section 3 Illustrates the current-mode controlled Boost converter example. Finally,

17.645797

18.35

12.044783

Semi−analytical simulation Scicos simulation

Simulation of Piecewise Hybrid Dynamical Systems in Matlab 5

11.096993

X(1) −2.0 −1.5 −1.0 −0.5 0.0 0.5 1.0 1.5

A general definition of HSs is presented here. This type of dynamical systems is characterized by the coexistence of two kinds of state vectors: continuous state vector *X*(*t*) of real values,

*<sup>X</sup>*˙ (*t*) = *<sup>F</sup>* (*X*(*t*), *<sup>E</sup>*(*t*)), *<sup>F</sup>* : H → **<sup>R</sup>***<sup>n</sup>*

<sup>H</sup> <sup>=</sup> **<sup>R</sup>***<sup>n</sup>* × M *is called the hybrid state space. X*(*t*) <sup>∈</sup> **<sup>R</sup>***<sup>n</sup> is the continuous state vector of the HS at time instant t and E*(*t*) ∈ M :<sup>=</sup> {1, . . . , *<sup>N</sup>*} *is its discrete state. E*+(*t*) *denotes the updated discrete state right after time instant t. φ* : H→M *describes the discrete dynamic, it is usually modeled by Petri Nets. A transition from E*(*t*) = *i to E*+(*t*) = *j is valid when the state X reaches a switching set*

*. Such transitions are called state dependent events. A HS is called piece-wise affine if for*

*<sup>E</sup>*+(*t*) = *<sup>φ</sup>* (*X*(*t*), *<sup>E</sup>*(*t*)), *<sup>φ</sup>* : H→M (1)

a conclusion is stated in Section 4.

X(2)

−1.0

−1.5

**2. Main results**

*called* S*Ei*,*Ej*

**Figure 1.** Numerical simulation versus semi–analytical simulation.

and discrete state vector *E*(*t*) belonging to a countable discrete set M.

**Definition 1.** *A continuous-time, autonomous HS is a system of the form:*

*each E* ∈ M*, F*(*X*, *E*) *can be defined by F*(*X*, *E*) = *AEX* + *BE*, ∀ *X.*

**2.1. A HS** (*X*, *E*, *t*)**: general definition**

−0.5

0.0

0.5

1.0

1.5

2.0

Switched circuits behavior is mostly simulated by pure numerical methods where precision step is increased when the system is near a switching condition. Those numerical tools are widely used mainly because of their ease-of-use and their ability to simulate a wide range of circuits including nonlinear, time–variant, and non–autonomous systems.

Even if those simulators can reach the desired relative precision for a continuous trajectory, they can miss a switching condition and then diverge drastically from the trajectory as in figure 1. This could be annoying when one is interested by border collision bifurcations, or when local behavior is needed with a good accuracy. In those applications, an alternative is to write down analytical, or semi–analytical, trajectories and switching conditions to obtain a recurrence which is very accurate and fast to run. Building and adapting such *ad'hoc* simulators represent a lot of efforts and a risk of mistakes.

Generic and accurate simulators can be proposed if we are restricted to a certain class of systems. A simulation tool with no loss of events is proposed in [14] and [15] for PWAHSs defined on polytopes (finite regions that are bounded by hyperplanes). This class of PWA differential systems has been widely studied as a standard technique to approximate a range of nonlinear systems.

But closed polytopic partition of the state space does not allow simulation of most switching circuits where switching frontiers are mostly single affine constraints or time–dependent periodical events.

This chapter follows our previous study in [13], [12] and [11].

We focus on planar PWAHSs with such simple switching conditions which can model a family of switched planar circuits: bang–bang regulators, the Boost converter, the Charge-Pump Phase Locked Loop (CP-PLL), . . .

This class of systems has analytical trajectories that help to build fast algorithms with no loss of events. We propose a semi-analytical solver for hybrid systems which provides:


This chapter is organized as follows. Section 2 contains our main results. We describe the problem to be deal with, we introduce a general algorithm to solve planar HSs, we present the algorithm that detects events' occurrence and devote a subsection to our approach efficiency. Section 3 Illustrates the current-mode controlled Boost converter example. Finally, a conclusion is stated in Section 4.

**Figure 1.** Numerical simulation versus semi–analytical simulation.

#### **2. Main results**

2 Will-be-set-by-IN-TECH

The study of nonlinear dynamics of DC-DC converters started in 1984 by works of Brockett and Wood [4]. Since then, chaos and nonlinear dynamics in power electronics circuits have attracted different research groups around the world. Different nonlinear phenomena have been observed such as routes to chaos following the period doubling cascade [16], [5], [19], [20] and [23] or quasi-periodic phenomena [6], [7] and [8], besides border collision bifurcations

Switched circuits behavior is mostly simulated by pure numerical methods where precision step is increased when the system is near a switching condition. Those numerical tools are widely used mainly because of their ease-of-use and their ability to simulate a wide range of

Even if those simulators can reach the desired relative precision for a continuous trajectory, they can miss a switching condition and then diverge drastically from the trajectory as in figure 1. This could be annoying when one is interested by border collision bifurcations, or when local behavior is needed with a good accuracy. In those applications, an alternative is to write down analytical, or semi–analytical, trajectories and switching conditions to obtain a recurrence which is very accurate and fast to run. Building and adapting such *ad'hoc*

Generic and accurate simulators can be proposed if we are restricted to a certain class of systems. A simulation tool with no loss of events is proposed in [14] and [15] for PWAHSs defined on polytopes (finite regions that are bounded by hyperplanes). This class of PWA differential systems has been widely studied as a standard technique to approximate a range

But closed polytopic partition of the state space does not allow simulation of most switching circuits where switching frontiers are mostly single affine constraints or time–dependent

We focus on planar PWAHSs with such simple switching conditions which can model a family of switched planar circuits: bang–bang regulators, the Boost converter, the Charge-Pump

This class of systems has analytical trajectories that help to build fast algorithms with no loss

• A pure analytic method when all continuous parts of the system and switching conditions can be solved symbolically. This can be the case for the boost converter [3], [21], the second

• A mixed method using analytical trajectories and numerical computation of the switching instant when those solutions are transcendent. This has been used for the third order

This chapter is organized as follows. Section 2 contains our main results. We describe the problem to be deal with, we introduce a general algorithm to solve planar HSs, we

of events. We propose a semi-analytical solver for hybrid systems which provides:

• A pure numerical method when the system is nonlinear or non-planar;

CP-PLL [17]. It can also be the case for the Buck converter [10], [21], . . .

circuits including nonlinear, time–variant, and non–autonomous systems.

simulators represent a lot of efforts and a risk of mistakes.

This chapter follows our previous study in [13], [12] and [11].

order charge-pump phase locked loop [17], [22].

[23] and [2].

of nonlinear systems.

periodical events.

Phase Locked Loop (CP-PLL), . . .

#### **2.1. A HS** (*X*, *E*, *t*)**: general definition**

A general definition of HSs is presented here. This type of dynamical systems is characterized by the coexistence of two kinds of state vectors: continuous state vector *X*(*t*) of real values, and discrete state vector *E*(*t*) belonging to a countable discrete set M.

**Definition 1.** *A continuous-time, autonomous HS is a system of the form:*

$$\begin{array}{l} \dot{X}(t) = F\left(X(t), E(t)\right), \; F: \mathcal{H} \to \mathbb{R}^n\\ E^+(t) = \phi\left(X(t), E(t)\right), \; \phi: \mathcal{H} \to \mathcal{M} \end{array} \tag{1}$$

<sup>H</sup> <sup>=</sup> **<sup>R</sup>***<sup>n</sup>* × M *is called the hybrid state space. X*(*t*) <sup>∈</sup> **<sup>R</sup>***<sup>n</sup> is the continuous state vector of the HS at time instant t and E*(*t*) ∈ M :<sup>=</sup> {1, . . . , *<sup>N</sup>*} *is its discrete state. E*+(*t*) *denotes the updated discrete state right after time instant t. φ* : H→M *describes the discrete dynamic, it is usually modeled by Petri Nets. A transition from E*(*t*) = *i to E*+(*t*) = *j is valid when the state X reaches a switching set called* S*Ei*,*Ej . Such transitions are called state dependent events. A HS is called piece-wise affine if for each E* ∈ M*, F*(*X*, *E*) *can be defined by F*(*X*, *E*) = *AEX* + *BE*, ∀ *X.*

**Remark —** For non–autonomous HSs, the function *φ* can also depend on time *φ* (*X*, *E*, *t*) : **<sup>R</sup>***<sup>n</sup>* ×M× **<sup>R</sup>** → M. Then time dependent events can occur and validate a transition, such as periodic events.

#### **2.2. HSs class of interest**

We consider a two dimensional PWAHS (*X*(*t*) <sup>∈</sup> **<sup>R</sup>**2). *<sup>F</sup>* has then the affine piece-wise form, *<sup>F</sup>*(., .) is defined for each *<sup>E</sup>* ∈ M and *<sup>X</sup>* <sup>∈</sup> **<sup>R</sup>**<sup>2</sup> by *<sup>F</sup>*(*X*(*t*), *<sup>E</sup>*(*t*)) = *AE*(*t*)*<sup>X</sup>* <sup>+</sup> *BE*(*t*), where *AE*(*t*) <sup>∈</sup> **<sup>R</sup>**2<sup>×</sup> <sup>2</sup> and *BE*(*t*) <sup>∈</sup> **<sup>R</sup>**<sup>2</sup> are two matrices that depend on the discrete state *<sup>E</sup>*(*t*). Hence, a two dimensional PWAHS is a HS that take the form:

$$\begin{aligned} \dot{X}(t) &= A\_{E(t)} X + B\_{E(t)'}\\ E(t) &\in \mathcal{M} = \{1, 2, \dots, N\} \end{aligned} \tag{2}$$

**2.3. Event–driven simulation of PWAHSs**

assuming that the discrete state is constant *E*(*t* > *t*

Compute all events' dates *d*S*EiEj* and *d*P*EiEj*

, *d*P*EiEj* );

We consider the affine Cauchy problem in **R**2:

and *E*(*t*

**end**

+

**Data**: *tk*, *Xk*, *Ek*. **while** *t* < *tfin* **do**

*tk*<sup>+</sup><sup>1</sup> = min(*d*S*EiEj*

*li*. The function *fi*(*t*) = *N*�

*2.4.1. Analytical trajectories*

*where* **I** *is the identity matrix.*

infinite value.

problem can be formulated as follows:

Find the smallest *t*

e*tA* =

∞ ∑ *k*=0 *t <sup>k</sup> A<sup>k</sup>*

*Xk*<sup>+</sup><sup>1</sup> = *f* (*Xk*, *Ek*, *tk*<sup>+</sup>1); *Ek*<sup>+</sup><sup>1</sup> = *φ* (*Xk*, *Ek*, *tk*<sup>+</sup>1);

The simulation will compute the hybrid state from event to event. Knowing the states *X*(*tk*)

;

**Algorithm 1**: Algorithm computing the hybrid state at *tk*+1.

*X*˙ (*t*) = *AX*(*t*) + *B*, *t* > *t*<sup>0</sup>

*X*(*t*0) = *X*<sup>0</sup>

where *X*<sup>0</sup> is the initial value. We compute the smallest strictly positive time *t*∗

*<sup>i</sup>* such that *fi*(*t*<sup>∗</sup>

If *fi* does not have any strictly positive root or the last condition is not satisfied, *t*<sup>∗</sup>

**Definition 2.** *For any square matrix A of order n and t in* **R***, the exponential matrix* e*tA is defined by*

It is well known that the analytical trajectory *X*(*t*) of the initial value matrix differential equation (3) is given in terms of the exponential matrix and the variation of constants formula

*<sup>k</sup>*! <sup>=</sup> **<sup>I</sup>** <sup>+</sup> *tA* <sup>+</sup>

trajectory of *X*(*t*) intersects the fixed border *Bi* arriving from the part of the plan where *N*�

*<sup>i</sup>* ) = 0 <sup>∃</sup> *<sup>δ</sup>* <sup>&</sup>gt; 0, <sup>∀</sup>*<sup>t</sup>* <sup>∈</sup>

> *t* <sup>2</sup> *A*<sup>2</sup> 2! <sup>+</sup>

+ *<sup>k</sup>* ) = *E*(*t* *tk f* 

+

*<sup>i</sup>* .*X*(*t*) − *li* defines the guard condition for a border *Bi*. Thus, the

*t* ∗ *<sup>i</sup>* − *δ*, *t*<sup>∗</sup> *i* 

*t* <sup>3</sup> *A*<sup>3</sup> *X*(*tk*), *E*(*t*

Simulation of Piecewise Hybrid Dynamical Systems in Matlab 7

+ *k* ) 

*<sup>k</sup>* ). Then the following algorithm

*dt* + *X*(*tk* ),

(3)

*<sup>i</sup>* .*X* <

*<sup>i</sup>* so that the

*<sup>i</sup>* is given the

, *fi*(*t*) <sup>&</sup>lt; <sup>0</sup> (4)

3! <sup>+</sup> ... (5)

*<sup>k</sup>* ), one can compute the trajectory *<sup>X</sup>*(*<sup>t</sup>* <sup>&</sup>gt; *tk*) = *<sup>t</sup>*

runs the simulation determining the date at the next event as the smallest:

**2.4. Event detection occurrence: description and algorithm**

∗

We consider two kinds of events: state dependent events and periodic events.

The state dependent event transition S*EiEj* is defined by an affine state border of the form *N*� *ij*.*X* < *lij*. In this case an event can occur when the continuous state reaches the border of the set S*EiEj* = *<sup>X</sup>*(*t*) <sup>∈</sup> **<sup>R</sup>**<sup>2</sup> : *<sup>N</sup>*� *ij*.*<sup>X</sup>* <sup>≤</sup> *lij* .

Note that the set S*EiEj* is not polytopic in the sense that the domain is not the interior of a closed bounded polytope.

**Remark —** We consider, without loss of generality, the case where a transition occurs at time *d*S*EiEj* if and only if the state *X*(*d*S*EiEj* ) reaches a border of the set S*EiEj* from outside. Figure 2 defines a transition with the complimentary set <sup>S</sup>¯, which allows to detect the event in both directions. Both transitions can be met with the set *<sup>B</sup>* <sup>=</sup> S ∪ <sup>S</sup>¯. Periodic

**Figure 2.** Oriented polytopic state dependent transitions.

events are simply defined by time instants *t* = *d*P*EiEj* , where *d*P*EiEj* belongs to the set P*EiEj* = {*t* : *t* = *kT* + *ϕ*, *k* ∈ **N**} . *T* is the period and *ϕ* is the phase of such periodic events.

#### **2.3. Event–driven simulation of PWAHSs**

4 Will-be-set-by-IN-TECH

**Remark —** For non–autonomous HSs, the function *φ* can also depend on time *φ* (*X*, *E*, *t*) : **<sup>R</sup>***<sup>n</sup>* ×M× **<sup>R</sup>** → M. Then time dependent events can occur and validate a transition, such as

We consider a two dimensional PWAHS (*X*(*t*) <sup>∈</sup> **<sup>R</sup>**2). *<sup>F</sup>* has then the affine piece-wise form, *<sup>F</sup>*(., .) is defined for each *<sup>E</sup>* ∈ M and *<sup>X</sup>* <sup>∈</sup> **<sup>R</sup>**<sup>2</sup> by *<sup>F</sup>*(*X*(*t*), *<sup>E</sup>*(*t*)) = *AE*(*t*)*<sup>X</sup>* <sup>+</sup> *BE*(*t*), where *AE*(*t*) <sup>∈</sup> **<sup>R</sup>**2<sup>×</sup> <sup>2</sup> and *BE*(*t*) <sup>∈</sup> **<sup>R</sup>**<sup>2</sup> are two matrices that depend on the discrete state *<sup>E</sup>*(*t*). Hence,

*<sup>X</sup>*˙ (*t*) = *AE*(*t*)*<sup>X</sup>* + *BE*(*t*),

The state dependent event transition S*EiEj* is defined by an affine state border of the form

Note that the set S*EiEj* is not polytopic in the sense that the domain is not the interior of a

**Remark —** We consider, without loss of generality, the case where a transition occurs at

Figure 2 defines a transition with the complimentary set <sup>S</sup>¯, which allows to detect the event in both directions. Both transitions can be met with the set *<sup>B</sup>* <sup>=</sup> S ∪ <sup>S</sup>¯. Periodic

Xk+1

Xk

P*EiEj* = {*t* : *t* = *kT* + *ϕ*, *k* ∈ **N**} . *T* is the period and *ϕ* is the phase of such periodic events.

*ij*.*X* < *lij*. In this case an event can occur when the continuous state reaches the border of

We consider two kinds of events: state dependent events and periodic events.

 .

*ij*.*X* ≤ *lij*

X(t,Xk,E(tk+))

N2.X<=l2

N1.X<=l1

*<sup>E</sup>*(*t*) ∈ M <sup>=</sup> {1, 2, . . . , *<sup>N</sup>*} (2)

) reaches a border of the set S*EiEj* from outside.

, where *d*P*EiEj* belongs to the set

periodic events.

*N*�

the set S*EiEj* =

time *d*S*EiEj*

closed bounded polytope.

**2.2. HSs class of interest**

a two dimensional PWAHS is a HS that take the form:

*<sup>X</sup>*(*t*) <sup>∈</sup> **<sup>R</sup>**<sup>2</sup> : *<sup>N</sup>*�

if and only if the state *X*(*d*S*EiEj*

**Figure 2.** Oriented polytopic state dependent transitions.

events are simply defined by time instants *t* = *d*P*EiEj*

The simulation will compute the hybrid state from event to event. Knowing the states *X*(*tk*) and *E*(*t* + *<sup>k</sup>* ), one can compute the trajectory *<sup>X</sup>*(*<sup>t</sup>* <sup>&</sup>gt; *tk*) = *<sup>t</sup> tk f X*(*tk*), *E*(*t* + *k* ) *dt* + *X*(*tk* ), assuming that the discrete state is constant *E*(*t* > *t* + *<sup>k</sup>* ) = *E*(*t* + *<sup>k</sup>* ). Then the following algorithm runs the simulation determining the date at the next event as the smallest:

**Data**: *tk*, *Xk*, *Ek*. **while** *t* < *tfin* **do** Compute all events' dates *d*S*EiEj* and *d*P*EiEj* ; *tk*<sup>+</sup><sup>1</sup> = min(*d*S*EiEj* , *d*P*EiEj* ); *Xk*<sup>+</sup><sup>1</sup> = *f* (*Xk*, *Ek*, *tk*<sup>+</sup>1); *Ek*<sup>+</sup><sup>1</sup> = *φ* (*Xk*, *Ek*, *tk*<sup>+</sup>1); **end**

**Algorithm 1**: Algorithm computing the hybrid state at *tk*+1.

#### **2.4. Event detection occurrence: description and algorithm**

We consider the affine Cauchy problem in **R**2:

$$\begin{cases} \dot{X}(t) = AX(t) + B, \; t > t\_0 \\ X(t\_0) = X\_0 \end{cases} \tag{3}$$

where *X*<sup>0</sup> is the initial value. We compute the smallest strictly positive time *t*∗ *<sup>i</sup>* so that the trajectory of *X*(*t*) intersects the fixed border *Bi* arriving from the part of the plan where *N*� *<sup>i</sup>* .*X* < *li*. The function *fi*(*t*) = *N*� *<sup>i</sup>* .*X*(*t*) − *li* defines the guard condition for a border *Bi*. Thus, the problem can be formulated as follows:

$$\begin{aligned} \text{Find the smallest } t\_i^\* \text{ such that } \begin{cases} f\_i(t\_i^\*) = 0\\ \exists \ \delta > 0, \forall t \in \left] t\_i^\* - \delta, t\_i^\* \right[\ \ \_\prime f\_i(t) < 0 \end{cases} \end{aligned} \tag{4}$$

If *fi* does not have any strictly positive root or the last condition is not satisfied, *t*<sup>∗</sup> *<sup>i</sup>* is given the infinite value.

#### *2.4.1. Analytical trajectories*

**Definition 2.** *For any square matrix A of order n and t in* **R***, the exponential matrix* e*tA is defined by*

$$\mathbf{e}^{tA} = \sum\_{k=0}^{\infty} \frac{t^k A^k}{k!} = \mathbb{I} + tA + \frac{t^2 A^2}{2!} + \frac{t^3 A^3}{3!} + \dots \tag{5}$$

*where* **I** *is the identity matrix.*

It is well known that the analytical trajectory *X*(*t*) of the initial value matrix differential equation (3) is given in terms of the exponential matrix and the variation of constants formula

#### 6 Will-be-set-by-IN-TECH 8 MATLAB – A Fundamental Tool for Scienti c Computing and Engineering Applications – Volume 3 Simulation of Piecewise Hybrid Dynamical Systems in Matlab <sup>7</sup>

by the general integral form:

$$\mathbf{X}(t) = \mathbf{e}^{(t-t\_0)A} \mathbf{X}\_0 + \int\_{t\_0}^t \mathbf{e}^{(t-s)A} \mathbf{B} \mathbf{ds}.\tag{6}$$

roots of the derivative function *f* �

**Table 2.** Expressions of *f* �

the algorithm

*T* ← {0, *L*, ∞};

**end**

**end**

*t* <sup>∗</sup> ← ∞;

**Data**: *Ni*, *li*, *A*, *B*, *Xk*

Break;

**3. Matlab modelling**

the set *L* of ordered roots of *f* �

**Result**: construct the set *L*, compute *t*∗

**if** *f*(*T*(*i*)) < 0 & *f*(*T*(*i* + 1)) > 0 **then** *t*<sup>∗</sup> ← *solve*[*T*(*i*), *T*(*i* + 1)];

elements, to find a crossing point when it exists.

**3.1. Application: Current-mode controlled Boost converter**

a switching controller. The basic circuit is given in figure 3.

and a load resistance connected in parallel with the capacitor.

**for** *i* ← 1 **to** (*card*(*T*) − 1) **do**

*f* �

(*t*) expressed in Table 2. We can then compute analytically

Simulation of Piecewise Hybrid Dynamical Systems in Matlab 9

, with *a*<sup>4</sup> ∈ **C**<sup>∗</sup>

(*t*), those roots determines monotone intervals of *f*(*t*). The

(*t*) = ... *p*<sup>1</sup> ∈ **R**<sup>∗</sup> *p*<sup>1</sup> = 0

*p*<sup>2</sup> = 0 *a*<sup>2</sup> + *a*<sup>4</sup> *p*<sup>1</sup> e*p*1*<sup>t</sup> a*<sup>2</sup> + 2 *a*<sup>3</sup> *t*

*<sup>p</sup>*<sup>1</sup> <sup>=</sup> *<sup>p</sup>*<sup>2</sup> <sup>∈</sup> **<sup>R</sup>**<sup>∗</sup> (*a*<sup>4</sup> *<sup>p</sup>*<sup>1</sup> <sup>+</sup> *<sup>a</sup>*<sup>5</sup> <sup>+</sup> *<sup>a</sup>*<sup>5</sup> *<sup>p</sup>*<sup>1</sup> *<sup>t</sup>*) <sup>e</sup>*p*1*<sup>t</sup>*

(*t*) depending on the eigenvalues *p*<sup>1</sup> and *p*<sup>2</sup>

following algorithm is used to return the solution *t*∗ when it exists or the value ∞ if not.

**Remark —** When (*p*1, *p*2) ∈ **C**<sup>∗</sup> × **C**<sup>∗</sup> the set *L* is infinite: when the real part of *pi* is positive,

**Algorithm 2**: Algorithm computing *t*∗ when a solution is transcendent.

will end by finding a root. In the other case, the set *L* should be reduced to its three first

Our semi-analytical solver is composed of different main programs that define the studied affine system. First, we create the affine system given with a specifically chosen name. Then, we define the matrices *Ai* and *Bi*. After that, we give the switching borders with the sign of transitions and all necessary elements or we give the period if it is about a periodic event. Finally, we execute the simulation by specifying the initial state and the time of simulation.

A current-mode controlled Boost converter in open loop consists of two parts: a converter and

This converter is a second-order circuit comprising an inductor, a capacitor, a diode, a switch

*<sup>p</sup>*<sup>1</sup> <sup>=</sup> *<sup>p</sup>*<sup>2</sup> <sup>∈</sup> **<sup>C</sup>**<sup>∗</sup> *<sup>a</sup>*<sup>4</sup> *<sup>p</sup>*<sup>1</sup> <sup>e</sup>*p*1*<sup>t</sup>* <sup>+</sup> *<sup>a</sup>*<sup>4</sup> *<sup>p</sup>*1e*<sup>p</sup>*<sup>1</sup> *<sup>t</sup>*

*<sup>p</sup>*<sup>2</sup> <sup>∈</sup> **<sup>R</sup>**<sup>∗</sup> *<sup>a</sup>*<sup>4</sup> *<sup>p</sup>*<sup>1</sup> <sup>e</sup> *<sup>p</sup>*1*<sup>t</sup>* <sup>+</sup> *<sup>a</sup>*<sup>6</sup> *<sup>p</sup>*<sup>2</sup> <sup>e</sup> *<sup>p</sup>*2*<sup>t</sup> <sup>a</sup>*<sup>2</sup> <sup>+</sup> *<sup>a</sup>*<sup>6</sup> *<sup>p</sup>*<sup>2</sup> <sup>e</sup>*p*2*<sup>t</sup>*

When *A* is invertible, the above expression becomes linear:

$$X(t) + A^{-1}B = \mathbf{e}^{(t - t\_0)A} (X\_0 + A^{-1}B) \tag{7}$$

The analytical expression of the exponential matrix e*At* takes two forms depending on whether the eigenvalues *p*<sup>1</sup> and *p*<sup>2</sup> of the matrix *A* are equal or not:

If *p*<sup>1</sup> �= *p*2, then

$$\mathbf{e}^{tA} = \frac{(p\_1 \mathbb{I} - A^{\diamond})}{p\_1 - p\_2} \mathbf{e}^{p\_1 t} - \frac{(p\_2 \mathbb{I} - A^{\diamond})}{p\_1 - p\_2} \mathbf{e}^{p\_2 t} \tag{8}$$

If *p*<sup>1</sup> = *p*<sup>2</sup> = *p*, then

$$\mathbf{e}^{tA} = \left(\mathbb{I} + \left(p\mathbb{I} - A^{\circ}\right)t\right)\mathbf{e}^{pt} \tag{9}$$
  $\text{where } A = \begin{pmatrix} a\_{11} \ a\_{12} \\ a\_{21} \ a\_{22} \end{pmatrix}, A^{\circ} = \begin{pmatrix} a\_{22} & -a\_{12} \\ -a\_{21} \ a\_{11} \end{pmatrix} \text{ and } \mathbb{I} = \begin{pmatrix} 1 \ 0 \\ 0 \ 1 \end{pmatrix}.$ 

Using these expressions, we can determine the function *f*(*t*) of the problem (4) as follows:

$$f(t) = a\_1 + a\_2t + a\_3t^2 + (a\_4 + a\_5t) \,\mathrm{e}^{p\_1t} + a\_6 \,\mathrm{e}^{p\_2t}$$

where *ai* are real scalars.

Depending on the eigenvalues *p*<sup>1</sup> and *p*2, there are five cases that determine the values of the coefficients *ai* as shown in Table 1. **Remark —** Coefficients *ai* are real scalars that depend on


**Table 1.** Expressions of *f*(*t*) depending on the eigenvalues *p*<sup>1</sup> and *p*2.

the eigenvalues *p*<sup>1</sup> and *p*2, the initial point *Xk* and the border parameters are *Ni* and *li*.

In some cases, (*p*<sup>1</sup> = *p*<sup>2</sup> = 0, gray cell in Table 1) roots of *f*(*t*) can be found analytically and the problem is solved with machine precision.

In other cases, the solution can not be found with classical functions and then a numeric algorithm should be used. Using classical methods like Newton does not guaranty existence or convergence of the smallest positive root. To meet these conditions, let us use analytical

roots of the derivative function *f* � (*t*) expressed in Table 2. We can then compute analytically


**Table 2.** Expressions of *f* � (*t*) depending on the eigenvalues *p*<sup>1</sup> and *p*<sup>2</sup>

the set *L* of ordered roots of *f* � (*t*), those roots determines monotone intervals of *f*(*t*). The following algorithm is used to return the solution *t*∗ when it exists or the value ∞ if not.

**Remark —** When (*p*1, *p*2) ∈ **C**<sup>∗</sup> × **C**<sup>∗</sup> the set *L* is infinite: when the real part of *pi* is positive, the algorithm

```
Data: Ni, li, A, B, Xk
Result: construct the set L, compute t∗
T ← {0, L, ∞};
t
∗ ← ∞;
for i ← 1 to (card(T) − 1) do
    if f(T(i)) < 0 & f(T(i + 1)) > 0 then
        t∗ ← solve[T(i), T(i + 1)];
        Break;
    end
end
```
**Algorithm 2**: Algorithm computing *t*∗ when a solution is transcendent.

will end by finding a root. In the other case, the set *L* should be reduced to its three first elements, to find a crossing point when it exists.

## **3. Matlab modelling**

6 Will-be-set-by-IN-TECH

The analytical expression of the exponential matrix e*At* takes two forms depending on

 *t t*0

<sup>e</sup>*p*<sup>1</sup> *<sup>t</sup>* <sup>−</sup> (*p*2**<sup>I</sup>** <sup>−</sup> *<sup>A</sup>*◦) *p*<sup>1</sup> − *p*<sup>2</sup>

> 1 0 0 1 .

<sup>2</sup> + (*a*<sup>4</sup> + *a*<sup>5</sup> *t*) e*p*<sup>1</sup> *<sup>t</sup>* + *a*<sup>6</sup> e*p*<sup>2</sup> *<sup>t</sup>*

*X*(*t*) + *A*−1*B* = e(*t*−*t*0)*A*(*X*<sup>0</sup> + *A*−1*B*) (7)

<sup>e</sup>*tA* <sup>=</sup> (**<sup>I</sup>** <sup>+</sup> (*<sup>p</sup>* **<sup>I</sup>** <sup>−</sup> *<sup>A</sup>*◦) *<sup>t</sup>*) <sup>e</sup>*pt* (9)

2

, with *a*<sup>5</sup> = *a*<sup>4</sup> ∈ **C**<sup>∗</sup>

e(*t*−*s*)*AB*d*s*. (6)

e*p*2*<sup>t</sup>* (8)

*X*(*t*) = e(*t*−*t*0)*AX*<sup>0</sup> +

When *A* is invertible, the above expression becomes linear:

whether the eigenvalues *p*<sup>1</sup> and *p*<sup>2</sup> of the matrix *A* are equal or not:

<sup>e</sup>*tA* <sup>=</sup> (*p*1**<sup>I</sup>** <sup>−</sup> *<sup>A</sup>*◦) *p*<sup>1</sup> − *p*<sup>2</sup>

> *<sup>a</sup>*<sup>22</sup> <sup>−</sup>*a*<sup>12</sup> −*a*<sup>21</sup> *a*<sup>11</sup>

*f*(*t*) = *a*<sup>1</sup> + *a*<sup>2</sup> *t* + *a*<sup>3</sup> *t*

Using these expressions, we can determine the function *f*(*t*) of the problem (4) as follows:

Depending on the eigenvalues *p*<sup>1</sup> and *p*2, there are five cases that determine the values of the coefficients *ai* as shown in Table 1. **Remark —** Coefficients *ai* are real scalars that depend on

*f*(*t*) = *a*1+... *p*<sup>1</sup> ∈ **R**<sup>∗</sup> *p*<sup>1</sup> = 0

*<sup>p</sup>*<sup>1</sup> <sup>=</sup> *<sup>p</sup>*<sup>2</sup> <sup>∈</sup> **<sup>R</sup>**<sup>∗</sup> (*a*<sup>4</sup> <sup>+</sup> *<sup>a</sup>*<sup>5</sup> *<sup>t</sup>*)e*p*1*<sup>t</sup>*

the eigenvalues *p*<sup>1</sup> and *p*2, the initial point *Xk* and the border parameters are *Ni* and *li*.

In some cases, (*p*<sup>1</sup> = *p*<sup>2</sup> = 0, gray cell in Table 1) roots of *f*(*t*) can be found analytically and

In other cases, the solution can not be found with classical functions and then a numeric algorithm should be used. Using classical methods like Newton does not guaranty existence or convergence of the smallest positive root. To meet these conditions, let us use analytical

*<sup>p</sup>*<sup>1</sup> <sup>=</sup> *<sup>p</sup>*<sup>2</sup> <sup>∈</sup> **<sup>C</sup>**<sup>∗</sup> *<sup>a</sup>*4e*<sup>p</sup>*1*<sup>t</sup>* <sup>+</sup> *<sup>a</sup>*5e*<sup>p</sup>*2*<sup>t</sup>*

**Table 1.** Expressions of *f*(*t*) depending on the eigenvalues *p*<sup>1</sup> and *p*2.

the problem is solved with machine precision.

*<sup>p</sup>*<sup>2</sup> <sup>∈</sup> **<sup>R</sup>**<sup>∗</sup> *<sup>a</sup>*<sup>4</sup> <sup>e</sup> *<sup>p</sup>*1*<sup>t</sup>* <sup>+</sup> *<sup>a</sup>*<sup>6</sup> <sup>e</sup> *<sup>p</sup>*2*<sup>t</sup> <sup>a</sup>*<sup>2</sup> *<sup>t</sup>* <sup>+</sup> *<sup>a</sup>*<sup>6</sup> <sup>e</sup>*p*2*<sup>t</sup> p*<sup>2</sup> = 0 *a*<sup>2</sup> *t* + *a*<sup>4</sup> e*p*1*<sup>t</sup> a*<sup>2</sup> *t* + *a*<sup>3</sup> *t*

and **I** =

by the general integral form:

If *p*<sup>1</sup> �= *p*2, then

where *A* =

If *p*<sup>1</sup> = *p*<sup>2</sup> = *p*, then

 *a*<sup>11</sup> *a*<sup>12</sup> *a*<sup>21</sup> *a*<sup>22</sup>

where *ai* are real scalars.

 , *A*◦ =

> Our semi-analytical solver is composed of different main programs that define the studied affine system. First, we create the affine system given with a specifically chosen name. Then, we define the matrices *Ai* and *Bi*. After that, we give the switching borders with the sign of transitions and all necessary elements or we give the period if it is about a periodic event. Finally, we execute the simulation by specifying the initial state and the time of simulation.

### **3.1. Application: Current-mode controlled Boost converter**

A current-mode controlled Boost converter in open loop consists of two parts: a converter and a switching controller. The basic circuit is given in figure 3.

This converter is a second-order circuit comprising an inductor, a capacitor, a diode, a switch and a load resistance connected in parallel with the capacitor.

The general circuit operation is driven by the switching controller. It compares the inductor current *iL* with the reference current *Iref* and generates the on/off driving signal for the switch *S*. When *S* is on, the current builds up in the inductor.

When the inductor current *iL* reaches a reference value, the switch opens and the inductor current flows through the load and the diode. The switch is again closed at the arrival of the next clock pulse from a free running clock of period *T*.

The Boost converter controlled in current mode is modeled by an affine piece-wise hybrid system defined by the same sub-systems given in equation as follows:

$$\begin{aligned} S\_1: \dot{X}(t) &= \begin{bmatrix} \frac{-1}{RC} & 0\\ 0 & 0 \end{bmatrix} X(t) + \begin{bmatrix} 0\\ \frac{V\_{in}}{L} \end{bmatrix} \\\\ S\_2: \dot{X}(t) &= \begin{bmatrix} \frac{-1}{RC} & \frac{1}{C} \\ \frac{-1}{L} & 0 \end{bmatrix} X(t) + \begin{bmatrix} 0\\ \frac{V\_{in}}{L} \end{bmatrix}, \\\\ \mathbf{L} &= \begin{bmatrix} \mathbf{L} & \mathbf{D} \\ \mathbf{C} & \mathbf{C} \mathbf{D} \end{bmatrix} \end{aligned} \tag{10}$$

where *T* is the period of this periodic event. The different simulations are obtained using our

Simulation of Piecewise Hybrid Dynamical Systems in Matlab 11

Before performing any study of the observed bifurcations in this circuit, a numerical

%%%%----------------------------calcule\_balayage\_mod.m----------------

The following program calcule\_balayage\_mod.m is used to obtain the parametric plane:

%% File calculating points to display a figure of parametric plane

% You should specify the path of hybrid\_solver\_matlab

planar PWA solver.

% %

%

%

%

%

eps=1E-6; ordre\_max=15; x\_eps = 1e-5; Xmax=100; nb\_trans=500; nb\_inits=1;

a=ta(1); b=tb(1);

L=1.5e-3; T=100e-6; R=40; Vin=b; C=5e-6; Iref=a; AE1=[

%

clear all; close all;

% Example of boost % Save data in...

monfich=('data\_balais');

ta= 0.5:1.1/200:1.6; tb= 5:15/200:20;

%% Definition of BOOST %Iref changes and noted a %Vin changes and noted b

> -1/R/C 0 ; 0 0 ];

addpath('.\hybrid\_solver\_matlab');

simulation in the parametric plane is needed.

**Figure 3.** Boost converter controlled in current mode.

In the case of the Boost converter controlled in current mode, there are two types of events: A state event defined by a fixed border of the set S*E*1*E*<sup>2</sup> :

$$\mathcal{S}\_{E\_1 E\_2} = \left\{ X \in \mathbb{R}^2 : [0 \ 1] \ X < I\_{ref} \right\} \tag{11}$$

and another periodic event defined by the dates *t* = *d*P*EiEj* , where *d*P*EiEj* belongs to the set:

$$\mathcal{P}\_{E\_2E\_1} = \{ t \in \mathbb{R} : t = nT, \ n \in \mathbb{N} \} \tag{12}$$

where *T* is the period of this periodic event. The different simulations are obtained using our planar PWA solver.

8 Will-be-set-by-IN-TECH

The general circuit operation is driven by the switching controller. It compares the inductor current *iL* with the reference current *Iref* and generates the on/off driving signal for the switch

When the inductor current *iL* reaches a reference value, the switch opens and the inductor current flows through the load and the diode. The switch is again closed at the arrival of the

The Boost converter controlled in current mode is modeled by an affine piece-wise hybrid

*X*(*t*) +

*X*(*t*) +

D

 0 *Vin L*

> 0 *Vin L*

 ,

R C

P*E*2*E*<sup>1</sup> = {*t* ∈ **R** : *t* = *nT*, *n* ∈ **N**} (12)

, where *d*P*EiEj* belongs to the set:

(10)

(11)

vC

 −1 *RC* 0 0 0 

 −1 *RC* 1 *C* −1 *<sup>L</sup>* 0

S

S

<sup>+</sup> QR

Q'

In the case of the Boost converter controlled in current mode, there are two types of events:

*<sup>X</sup>* <sup>∈</sup> **<sup>R</sup>**<sup>2</sup> : [0 1] *<sup>X</sup>* <sup>&</sup>lt; *Iref*

*S*. When *S* is on, the current builds up in the inductor.

next clock pulse from a free running clock of period *T*.

iL L

iL

−

T

A state event defined by a fixed border of the set S*E*1*E*<sup>2</sup> :

and another periodic event defined by the dates *t* = *d*P*EiEj*

S*E*1*E*<sup>2</sup> =

**Figure 3.** Boost converter controlled in current mode.

Horloge

Iref

+

Vin

system defined by the same sub-systems given in equation as follows:

*S*<sup>1</sup> : *X*˙ (*t*) =

*S*<sup>2</sup> : *X*˙ (*t*) =

Before performing any study of the observed bifurcations in this circuit, a numerical simulation in the parametric plane is needed.

The following program calcule\_balayage\_mod.m is used to obtain the parametric plane:

```
%%%%----------------------------calcule_balayage_mod.m----------------
%
%
clear all;
close all;
%
%% File calculating points to display a figure of parametric plane
% Example of boost
% Save data in...
%
monfich=('data_balais');
% You should specify the path of hybrid_solver_matlab
%
addpath('.\hybrid_solver_matlab');
%
eps=1E-6;
ordre_max=15;
x_eps = 1e-5;
Xmax=100;
nb_trans=500;
nb_inits=1;
ta= 0.5:1.1/200:1.6;
tb= 5:15/200:20;
a=ta(1);
b=tb(1);
%% Definition of BOOST
%Iref changes and noted a
%Vin changes and noted b
%
L=1.5e-3;
T=100e-6;
R=40;
Vin=b;
C=5e-6;
Iref=a;
AE1=[
    -1/R/C 0 ;
    0 0
    ];
```
10 Will-be-set-by-IN-TECH 12 MATLAB – A Fundamental Tool for Scienti c Computing and Engineering Applications – Volume 3 Simulation of Piecewise Hybrid Dynamical Systems in Matlab <sup>11</sup>

break;

or=0; break;

if (or==ordre\_max)

if (X.E ~= 1)

it=1; X0 = X; tt0=X.t;

[X]=recu(H,X);

tt0=X.t;

else

end

end

end or;

end

end

temps=toc save(monfich) affiche\_balayage

end

or=it; break;

if (max(abs(X.Xc))>Xmax)

%% check if we have periodic event state E=1

[X]=recu(H,X); %1->2;

ii=1;

end [X]=recu(H,X);

it = it + 1;

ordres(ib,ia)=max(or,ordres(ib,ia));

while (it<ordre\_max) & (or == ordre\_max)

Simulation of Piecewise Hybrid Dynamical Systems in Matlab 13

while (ii\*T<tt) it=it+1; ii=ii+1;

if (max(abs(X.Xc-X0.Xc))<x\_eps)

fprintf('About %2.1f %% done, still about %5.0f secondes to be... %3.0f minutes\n',ia/na\*100,toc/ia\*(na-ia),toc/ia\*(na-ia)/60)

After calculating the necessary points of the parametric plane saved in the file named

dat\_balais, the next program affiche\_balayage.m plots the figure given in Fig.4

tt=X.t-tt0;

end

end

end

end

```
AE2=[
    -1/R/C 1/C ;
    -1/L 0
    ];
B1=[0;
    b/L];
B2=B1;
N1 = [0 1];
S1 = '<';
[p1,p2]=racines(AE1);
H=create_hybrid_system('affine');
H=add_state(H,1,'On',AE1,B1);
H=add_state(H,2,'Off',AE2,B2);
H=add_event(H,1,'Iref',N1,a,S1);
H=add_periodic_event(H,1,'Clock',T,0);
H=add_transition(H,1,2,1);
H=add_periodic_transition(H,2,1,1);
%
%% initial state
Xi.t=T/1000;
Xi.E=1;
Xi.Xc=[16.5;0.47];
mape=colormap(ma_color);
na=length(ta)
tic
for ia = 1 : na
    a=ta(ia);
    for ib= 1 : length(tb)
        b=tb(ib);
        %% update of the equation with a new a
        % here only Iref that changes and modifies a border
        B1=[0;
            b/L];
        H=add_state(H,1,'On',AE1,B1);
        H=add_state(H,2,'Off',AE2,B2);
        H=add_event(H,1,'Iref',N1,a,S1);
        H=update_transition(H,1,2,1);
        ordres(ib,ia)=-2;
        for init = 1:nb_inits
            X=Xi;
             or = ordre_max;
            for i = 1 : nb_trans
                [X]=recu(H,X); %1->2;
                [X]=recu(H,X); %2->1;
                if (X.t==Inf)
                    or=-1;
```

```
break;
                 end
                 if (max(abs(X.Xc))>Xmax)
                     or=0;
                     break;
                 end
            end
            if (or==ordre_max)
      %% check if we have periodic event state E=1
               if (X.E ~= 1)
                  [X]=recu(H,X);
               end
                 it=1;
                 X0 = X;
                 tt0=X.t;
                 while (it<ordre_max) & (or == ordre_max)
                     [X]=recu(H,X); %1->2;
                              tt=X.t-tt0;
                              ii=1;
                              while (ii*T<tt)
                                  it=it+1;
                                  ii=ii+1;
                              end
                     [X]=recu(H,X);
                     tt0=X.t;
                     if (max(abs(X.Xc-X0.Xc))<x_eps)
                         or=it;
                         break;
                     else
                         it = it + 1;
                     end
                 end
            end
            or;
            ordres(ib,ia)=max(or,ordres(ib,ia));
        end
    end
    fprintf('About %2.1f %% done, still about %5.0f secondes to be...
     %3.0f minutes\n',ia/na*100,toc/ia*(na-ia),toc/ia*(na-ia)/60)
end
temps=toc
save(monfich)
affiche_balayage
```
10 Will-be-set-by-IN-TECH

AE2=[

B2=B1; N1 = [0 1]; S1 = '<';

%

tic

]; B1=[0;

b/L];


[p1,p2]=racines(AE1);

H=add\_transition(H,1,2,1);

%% initial state Xi.t=T/1000; Xi.E=1;

Xi.Xc=[16.5;0.47];

na=length(ta)

for ia = 1 : na a=ta(ia);

mape=colormap(ma\_color);

B1=[0;

for ib= 1 : length(tb) b=tb(ib);

b/L];

ordres(ib,ia)=-2; for init = 1:nb\_inits

X=Xi;

%% update of the equation with a new a

[X]=recu(H,X); %1->2; [X]=recu(H,X); %2->1;

H=add\_state(H,1,'On',AE1,B1); H=add\_state(H,2,'Off',AE2,B2); H=add\_event(H,1,'Iref',N1,a,S1); H=update\_transition(H,1,2,1);

> or = ordre\_max; for i = 1 : nb\_trans

> > if (X.t==Inf) or=-1;

% here only Iref that changes and modifies a border

H=create\_hybrid\_system('affine'); H=add\_state(H,1,'On',AE1,B1); H=add\_state(H,2,'Off',AE2,B2); H=add\_event(H,1,'Iref',N1,a,S1);

H=add\_periodic\_event(H,1,'Clock',T,0);

H=add\_periodic\_transition(H,2,1,1);

After calculating the necessary points of the parametric plane saved in the file named dat\_balais, the next program affiche\_balayage.m plots the figure given in Fig.4

12 Will-be-set-by-IN-TECH 14 MATLAB – A Fundamental Tool for Scienti c Computing and Engineering Applications – Volume 3 Simulation of Piecewise Hybrid Dynamical Systems in Matlab <sup>13</sup>

The figure 4 of the parametric plane allows to emphasize the parameters values for which

Simulation of Piecewise Hybrid Dynamical Systems in Matlab 15

Figure 5 shows a bifurcation diagram (Feigenbaum type) in the plane (*Iref* , *iL*). However,

there exists at least one attractor (fixed point, cycle of order *k*, strange attractor).

To draw these two figures we use programs: calcule\_figuier.m and affiche\_figuier.m

%%%%----------------------------calcule\_figuier.m----------------

%% Code that calculates then displays the points of a bifurcation

nb\_trans=400;%400 %Number of iterations to pass the transient phase ordre\_max=100; %100% nombre de points affichés après le transitoire ta= 0.5:0.0025:1.6; % values of the parameter a to be calculated

% points (x or y, index a, the ordre\_max of the last trajectory

figure 6 shows the bifurcation diagram in the space plane (*I*, *iL*, *vC*).

%

%

%

%

points) a=ta(1);

L=1.5e-3; T=100e-6; R=40; Vin=10; C=5e-6; Iref=a; AE1=[

]; AE2=[

];



%

tree

clear all; close all;

% Boost converter example % Save the data in ... monfich=('data\_points');

addpath('../hybrid\_solver\_matlab');

eps=1E-6; % precision of the solver

%ta= 1.22:0.001:1.4;%0.5:0.001:1.6; points=zeros(2,length(ta),ordre\_max);

%% Definition of the Boost converter

%Iref is a variable denoted a

```
%%%%----------------------------affiche_balayage_mod.m-------------
% Used in general after calcule_balayage
%% Charge the saved 2-D bifurcation scan
%if the file was not executed
if (exist('ordres')==1)
    disp('use the points matrix of the workspace');
elseif (exist('data_balais.mat')==2)
    disp('charge the points that are in data_balais.mat');
    load data_balais.mat
else
    disp('There are no points or files of points: try ordres.mat...
           insha ALLAH! It may be long...')
    load ordres.mat
end
%
%% Display the bifurcation scan diagram
da=(ta(2)-ta(1))/2;
db=(tb(2)-tb(1))/2;
colormap(mape)
for ia = 1 : length(ta)
    a=ta(ia);
    for ib= 1 : length(tb)
        b=tb(ib);
        if (ordres(ib,ia)<0)
            %plot(a,b,'.w');
            fill([a-da a-da a+da a+da],[b-db b+db b+db b-db],...
            'w','EdgeColor','none')
        elseif (ordres(ib,ia)==0)
            %plot(a,b,'+w');
            fill([a-da a-da a+da a+da],[b-db b+db b+db b-db],...
            'w','EdgeColor','none')
        else
            %plot(a,b,'s','color',mape(ordres(ib,ia),:),...
            'MarkerFaceColor',mape(ordres(ib,ia),:),'MarkerSize');
            fill([a-da a-da a+da a+da],[b-db b+db b+db b-db],...
            mape(ordres(ib,ia),:),'EdgeColor','none')
        end
       hold on
    end
end
colormap(mape)
colorbar %(mape)
colorbar('YTickLabel',...
        {'O1','O2','O3','O4','O5','O6','O7',...
         'O8','O9','10','11','12','13','14','O+'})
```
The figure 4 of the parametric plane allows to emphasize the parameters values for which there exists at least one attractor (fixed point, cycle of order *k*, strange attractor).

Figure 5 shows a bifurcation diagram (Feigenbaum type) in the plane (*Iref* , *iL*). However, figure 6 shows the bifurcation diagram in the space plane (*I*, *iL*, *vC*).

To draw these two figures we use programs: calcule\_figuier.m and affiche\_figuier.m

12 Will-be-set-by-IN-TECH

%%%%----------------------------affiche\_balayage\_mod.m-------------

disp('use the points matrix of the workspace');

insha ALLAH! It may be long...')

disp('charge the points that are in data\_balais.mat');

disp('There are no points or files of points: try ordres.mat...

fill([a-da a-da a+da a+da],[b-db b+db b+db b-db],...

fill([a-da a-da a+da a+da],[b-db b+db b+db b-db],...

'MarkerFaceColor',mape(ordres(ib,ia),:),'MarkerSize'); fill([a-da a-da a+da a+da],[b-db b+db b+db b-db],...

%plot(a,b,'s','color',mape(ordres(ib,ia),:),...

mape(ordres(ib,ia),:),'EdgeColor','none')

{'O1','O2','O3','O4','O5','O6','O7',... 'O8','O9','10','11','12','13','14','O+'})

% Used in general after calcule\_balayage

%% Charge the saved 2-D bifurcation scan

elseif (exist('data\_balais.mat')==2)

%% Display the bifurcation scan diagram

%if the file was not executed

load data\_balais.mat

if (exist('ordres')==1)

load ordres.mat

for ia = 1 : length(ta)

for ib= 1 : length(tb) b=tb(ib);

> if (ordres(ib,ia)<0) %plot(a,b,'.w');

> > 'w','EdgeColor','none')

'w','EdgeColor','none')

elseif (ordres(ib,ia)==0) %plot(a,b,'+w');

da=(ta(2)-ta(1))/2; db=(tb(2)-tb(1))/2; colormap(mape)

a=ta(ia);

else

end hold on

colorbar('YTickLabel',...

end

colormap(mape) colorbar %(mape)

end

else

end %

```
%%%%----------------------------calcule_figuier.m----------------
%
clear all;
close all;
%
%% Code that calculates then displays the points of a bifurcation
tree
% Boost converter example
% Save the data in ...
monfich=('data_points');
%
addpath('../hybrid_solver_matlab');
%
eps=1E-6; % precision of the solver
nb_trans=400;%400 %Number of iterations to pass the transient phase
ordre_max=100; %100% nombre de points affichés après le transitoire
ta= 0.5:0.0025:1.6; % values of the parameter a to be calculated
%ta= 1.22:0.001:1.4;%0.5:0.001:1.6;
points=zeros(2,length(ta),ordre_max);
% points (x or y, index a, the ordre_max of the last trajectory
points)
a=ta(1);
%
%% Definition of the Boost converter
%Iref is a variable denoted a
L=1.5e-3;
T=100e-6;
R=40;
Vin=10;
C=5e-6;
Iref=a;
AE1=[
    -1/R/C 0 ;
       0 0
    ];
AE2=[
    -1/R/C 1/C ;
    -1/L 0
    ];
```
Xin=X.Xc; tt0=X.t; it=1;

while (it<(ordre\_max+1)) [X]=recu(H,X);%1->2

tt=X.t-tt0;

while (ii\*T<tt)

it=it+1; ii=ii+1;

[X]=recu(H,X);%2->1

(na-ia)/60)

%scan the values of 'a'

if (exist('points')==1)

load data\_points

points(:,ia,it)=X.Xc;

ii=1;

end

Xin=X.Xc; tt0=X.t; it=it+1;

%

end %

end

%% cc vi

seconds...

temps\_ecoule=toc save(monfich)

affiche\_figuier

2-> 1 %

%% memorize ordremax points issued from the periodic transition

Simulation of Piecewise Hybrid Dynamical Systems in Matlab 17

points(:,ia,it)=traj\_ni(AE1,B1,p1,p2,Xin,ii\*T);

fprintf('Approximately %2.1f %% are done, yet approximately %5.0f

%%%%----------------------------affiche\_figuier.m----------------

%% Charges the file containing the saved points

disp('use the points matrix of the workspace');

disp('charge the points that are in data\_points.mat');

% if the file "figuier" is not executed

elseif (exist('data\_points.mat')==2)

that is %3.0f minutes\n',ia/na\*100,toc/ia\*(na-ia),toc/ia\*

Xin=traj\_ni(AE1,B1,p1,p2,Xin,ii\*T);

```
B1=[0;
    Vin/L];
B2=B1;
N1 = [0 1];
S1 = '<';
[p1,p2]=racines(AE1)
H=create_hybrid_system('affine');
H=add_state(H,1,'On',AE1,B1);
H=add_state(H,2,'Off',AE2,B2);
H=add_event(H,1,'Iref',N1,a,S1);
H=add_periodic_event(H,2,'Clock',T,0);
H=add_transition(H,1,2,1);
H=add_periodic_transition(H,2,1,2);
%
%% Initial condition
X0.t=T/1000;
X0.E=1;
X0.Xc=[16.4549;0.4648];
%
%% Vary a and memorize the points for the bifurcation tree
na=length(ta);
tic;
for ia = 1 : na
    a=ta(ia);
    if a==1.3
    end
    %% Update of the equation with a new a
     vi=X.Xc;
     cc=ia;
    % Here only Iref varies and the corresponding border is then
    modified
    H=add_event(H,1,'Iref',N1,a,S1);
    H=update_transition(H,1,2,1);
    X=X0;
    %
    %% transient zone
    for i = 1: nb_trans
        [X]=recu(H,X);
          [X]=recu(H,X);
    end
    %% Assure that we are on a periodic event, state E=1
    %
     if (X.E ~= 1)
         [X]=recu(H,X);
     end
     %
```

```
Xin=X.Xc;
       tt0=X.t;
       it=1;
    %% memorize ordremax points issued from the periodic transition
    2-> 1
    %
     while (it<(ordre_max+1))
        [X]=recu(H,X);%1->2
        tt=X.t-tt0;
        ii=1;
        while (ii*T<tt)
            points(:,ia,it)=traj_ni(AE1,B1,p1,p2,Xin,ii*T);
            Xin=traj_ni(AE1,B1,p1,p2,Xin,ii*T);
            it=it+1;
            ii=ii+1;
        end
      %
      [X]=recu(H,X);%2->1
      points(:,ia,it)=X.Xc;
      Xin=X.Xc;
      tt0=X.t;
      it=it+1;
    end
   %
  fprintf('Approximately %2.1f %% are done, yet approximately %5.0f
  seconds...
          that is %3.0f minutes\n',ia/na*100,toc/ia*(na-ia),toc/ia*
          (na-ia)/60)
end
%scan the values of 'a'
temps_ecoule=toc
save(monfich)
%%
cc
vi
affiche_figuier
%%%%----------------------------affiche_figuier.m----------------
%% Charges the file containing the saved points
% if the file "figuier" is not executed
if (exist('points')==1)
   disp('use the points matrix of the workspace');
elseif (exist('data_points.mat')==2)
   disp('charge the points that are in data_points.mat');
   load data_points
```
14 Will-be-set-by-IN-TECH

B1=[0;

B2=B1; N1 = [0 1]; S1 = '<';

%

%

tic;

Vin/L];

[p1,p2]=racines(AE1)

%% Initial condition

X0.Xc=[16.4549;0.4648];

X0.t=T/1000; X0.E=1;

na=length(ta);

for ia = 1 : na a=ta(ia); if a==1.3

> vi=X.Xc; cc=ia;

modified

%% transient zone for i = 1: nb\_trans [X]=recu(H,X); [X]=recu(H,X);

if (X.E ~= 1)

[X]=recu(H,X);

X=X0; %

end

end %

%

end

H=add\_transition(H,1,2,1);

H=create\_hybrid\_system('affine'); H=add\_state(H,1,'On',AE1,B1); H=add\_state(H,2,'Off',AE2,B2); H=add\_event(H,1,'Iref',N1,a,S1);

H=add\_periodic\_event(H,2,'Clock',T,0);

%% Vary a and memorize the points for the bifurcation tree

% Here only Iref varies and the corresponding border is then

%% Assure that we are on a periodic event, state E=1

%% Update of the equation with a new a

H=add\_event(H,1,'Iref',N1,a,S1); H=update\_transition(H,1,2,1);

H=add\_periodic\_transition(H,2,1,2);

for ia=1:length(ta)

xlabel('Iref(A)'); ylabel('iL(A)'); zlabel('vC(V)');

associated with one color.

draw bifurcation curves.

%%

%

%

%

% Warning

end

end

for io=1:ordre\_max

plot3(ta(ia),points(2,ia,io),points(1,ia,io));

Simulation of Piecewise Hybrid Dynamical Systems in Matlab 19

In these two figures, the voltage *Vin* is fixed to 10V and the current *Iref* varies in the interval [0.5, 1.6]. We observe a period cascade doubling leading to a chaotic regime, interrupted by a border collision bifurcation at *Iref* = 1.23A (see figure 7). In this figure, a distinction is given between the attractors of attractive cycle type of the order 1 to 14. Each cycle of order *k* is

For example, the blue area O1 represents the parameters' values for which there exists an attractive fixed point (fundamental periodic regime). The red area O2 represents the existence of an attractive cycle of order 2. The yellow area O4 represents the existence of an attractive cycle of order 4 and so on until getting the cycles O14 of order *k* = 14. The black area O+ corresponds to parameters values (*Iref* , *Vin*) for which there exist cycles of order *k* ≥ 15 or other types of attractors. In this last area, a chaotic phenomenon could be observed. This bi-dimensional diagram shows some bifurcation curves. In fact, for the rectangle defined by the interval of parameter *Vin* ∈ [7, 15] and the parameter *Iref* ∈ [0.5, 1.6], we observe an area of blue color (existence of attractive fixed point) followed by an area of red color (existence of cycle of order 2), an area of yellow color (existence of cycle of order 4) and another area of black color (existence of cycle of order *k* ≥ 15 or another attractor type); this succession of

This representation of the parametric plane is not enough to establish a bifurcation structure of the hybrid model of the Boost converter, but it is useful for the initialization of programs to

The simulation results (temporal domain and voltage-current plane (*vC*, *iL*)) are obtained using the planar PWH solver in the case of the Boost converter controlled in current mode for periods: 1*T* (figure 8) for *Iref* = 0.7A , 2*T* (figure 9) for *Iref* = 1A, 4*T* (figure 10) for *Iref* = 1.3A and the chaotic regime (figure 11) for *Iref* = 1.5A. For these plots we used the code below by choosing the bifurcation parameter *Iref* corresponding to each period case.

zones corresponds to the existence of period doubling cascade.

% Define an affine system ii at random and simulate it

% In the affine case '=' is not supported yet !

% by use of our Matlab toolbox solver

%detection of errors manually

addpath('.\hybrid\_solver');

```
else
 disp('There are no points or files of points: try points.mat insh
 ALLAH!')
 load points
end
%% bifurcation tree depending on the dimensions x then y
% % for dim=1:2
% %
% % plot(points(dim,1,1));
% % hold on;
% % for ia=1:length(ta)
% % for io=2:ordre_max
% % plot(ta(ia),points(dim,ia,io));
% % end
% % end
% % xlabel('a');
% % figure
% % end
% figure
%
% plot(points(1,1,1));
% hold on;
% for ia=1:length(ta)
% for io=2:ordre_max
% plot(ta(ia),points(1,ia,io));
% end
% end
% xlabel('Vin(V)');
% ylabel('vC(V)');
   plot(points(2,1,1));
   hold on;
   for ia=1:length(ta)
       for io=2:ordre_max
           plot(ta(ia),points(2,ia,io));
       end
   end
   xlabel('Iref(A)');
   ylabel('iL(A)');
%% bifurcation tree in 3D
% z = variable parameter denoted a
% x the dimension x
% y the dimension y of the point
 plot3(points(1,1,1),points(2,1,1),ta(1));
 plot3(ta(1),points(2,1,1),points(1,1,1));
 hold on;
```

```
for ia=1:length(ta)
    for io=1:ordre_max
        plot3(ta(ia),points(2,ia,io),points(1,ia,io));
    end
end
xlabel('Iref(A)');
ylabel('iL(A)');
zlabel('vC(V)');
```
16 Will-be-set-by-IN-TECH

disp('There are no points or files of points: try points.mat insh

%% bifurcation tree depending on the dimensions x then y

% % plot(ta(ia),points(dim,ia,io));

% plot(ta(ia),points(1,ia,io));

plot(ta(ia),points(2,ia,io));

else

end

% %

ALLAH!') load points

% % for dim=1:2

% % hold on;

% % end % % end

% % figure

% hold on;

% end % end

hold on;

end

% x the dimension x

end

hold on;

% % end % figure

%

% % xlabel('a');

% % plot(points(dim,1,1));

% % for ia=1:length(ta) % % for io=2:ordre\_max

% plot(points(1,1,1));

% for ia=1:length(ta) % for io=2:ordre\_max

% xlabel('Vin(V)'); % ylabel('vC(V)');

plot(points(2,1,1));

for ia=1:length(ta)

xlabel('Iref(A)'); ylabel('iL(A)'); %% bifurcation tree in 3D

for io=2:ordre\_max

% z = variable parameter denoted a

plot3(points(1,1,1),points(2,1,1),ta(1)); plot3(ta(1),points(2,1,1),points(1,1,1));

% y the dimension y of the point

In these two figures, the voltage *Vin* is fixed to 10V and the current *Iref* varies in the interval [0.5, 1.6]. We observe a period cascade doubling leading to a chaotic regime, interrupted by a border collision bifurcation at *Iref* = 1.23A (see figure 7). In this figure, a distinction is given between the attractors of attractive cycle type of the order 1 to 14. Each cycle of order *k* is associated with one color.

For example, the blue area O1 represents the parameters' values for which there exists an attractive fixed point (fundamental periodic regime). The red area O2 represents the existence of an attractive cycle of order 2. The yellow area O4 represents the existence of an attractive cycle of order 4 and so on until getting the cycles O14 of order *k* = 14. The black area O+ corresponds to parameters values (*Iref* , *Vin*) for which there exist cycles of order *k* ≥ 15 or other types of attractors. In this last area, a chaotic phenomenon could be observed. This bi-dimensional diagram shows some bifurcation curves. In fact, for the rectangle defined by the interval of parameter *Vin* ∈ [7, 15] and the parameter *Iref* ∈ [0.5, 1.6], we observe an area of blue color (existence of attractive fixed point) followed by an area of red color (existence of cycle of order 2), an area of yellow color (existence of cycle of order 4) and another area of black color (existence of cycle of order *k* ≥ 15 or another attractor type); this succession of zones corresponds to the existence of period doubling cascade.

This representation of the parametric plane is not enough to establish a bifurcation structure of the hybrid model of the Boost converter, but it is useful for the initialization of programs to draw bifurcation curves.

The simulation results (temporal domain and voltage-current plane (*vC*, *iL*)) are obtained using the planar PWH solver in the case of the Boost converter controlled in current mode for periods: 1*T* (figure 8) for *Iref* = 0.7A , 2*T* (figure 9) for *Iref* = 1A, 4*T* (figure 10) for *Iref* = 1.3A and the chaotic regime (figure 11) for *Iref* = 1.5A. For these plots we used the code below by choosing the bifurcation parameter *Iref* corresponding to each period case.

```
%%
% Define an affine system ii at random and simulate it
% by use of our Matlab toolbox solver
%
%detection of errors manually
%
% Warning
% In the affine case '=' is not supported yet !
addpath('.\hybrid_solver');
%
```

```
L=1.5E-3;
T=100E-6;
R=40;
E=10;
C=5E-6;
%Bifurcation parameter
Iref=0.7; %1;%1.3;%1.5;
% X = [vc ; iL]
%
%System 1 On
A1=[ -1/R/C, 0 ;
     0 0 ];
B1=[0 ; E/L];
N1 = [0 1];
Lim1= Iref;
% System 2 Off
A2=[ -1/R/C, 1/C ;
     -1/L 0 ];
B2=[0 ; E/L];
%
%% Initial condition
X0.t=0;
X0.E=1;
X0.Xc=[16.4549;0.4648];%[13.6097 ;0.3435];
%
%% Boost converter
clear H;
H=create_hybrid_system('Boost Converter');
H=add_state(H,1,'On',A1,B1);
H=add_state(H,2,'Off',A2,B2);
%
H=add_event(H,1,'Iref',N1,Lim1,'<');
H=add_transition(H,1,2,1);
H=add_periodic_event(H,2,'Clock',T,0);
H=add_periodic_transition(H,2,1,2);
Han=H;
%
Xan = hsim(Han,X0,4*T);
%
[XcAn,EAn,tAn]=split_state(Xan);
%
 trajplane(Xan,Han)
     figure;
    subplot(211);
    trajplot(Xan,Han,1);
    subplot(212);
    trajplot(Xan,Han,2);
```
0.6 0.8 <sup>1</sup> 1.2 1.4 1.6 <sup>5</sup>

O1 O2 O3 O4 O5 O6 O7 O8 O9 O10 O11 O12 O13 O14 O+

Simulation of Piecewise Hybrid Dynamical Systems in Matlab 21

Iref(A)

**Figure 4.** Parametric diagram of the Boost converter in the plane (*Iref* , *Vin*) for *Iref* ∈ [0.5, 1.6]A and

0.5 0.6 0.8 <sup>1</sup> 1.2 1.4 1.6 0.2

Iref(A)

**Figure 5.** Bifurcation diagram of the Boost in the plane (*Iref* , *iL*) for *Iref* ∈ [0.5, 1.6]A and *Vin* = 10V fixed.

10

0.4

0.6

0.8

iL(A)

1

1.2

1.4

1.6

15

Vin(V)

*Vin* ∈ [5, 20]V.

20

18 Will-be-set-by-IN-TECH

L=1.5E-3; T=100E-6; R=40; E=10; C=5E-6;

%

%

%

%

%

%

Han=H; %

X0.t=0; X0.E=1;

clear H;

%Bifurcation parameter Iref=0.7; %1;%1.3;%1.5;

0 0 ];


H=add\_state(H,1,'On',A1,B1); H=add\_state(H,2,'Off',A2,B2);

H=add\_transition(H,1,2,1);

Xan = hsim(Han,X0,4\*T);

trajplane(Xan,Han) figure; subplot(211);

subplot(212);

X0.Xc=[16.4549;0.4648];%[13.6097 ;0.3435];

H=create\_hybrid\_system('Boost Converter');

H=add\_event(H,1,'Iref',N1,Lim1,'<');

[XcAn,EAn,tAn]=split\_state(Xan);

trajplot(Xan,Han,1);

trajplot(Xan,Han,2);

H=add\_periodic\_event(H,2,'Clock',T,0); H=add\_periodic\_transition(H,2,1,2);

%% Initial condition

%% Boost converter

% X = [vc ; iL]

%System 1 On A1=[ -1/R/C, 0 ;

B1=[0 ; E/L]; N1 = [0 1]; Lim1= Iref;

% System 2 Off A2=[ -1/R/C, 1/C ;

B2=[0 ; E/L];

**Figure 4.** Parametric diagram of the Boost converter in the plane (*Iref* , *Vin*) for *Iref* ∈ [0.5, 1.6]A and *Vin* ∈ [5, 20]V.

**Figure 5.** Bifurcation diagram of the Boost in the plane (*Iref* , *iL*) for *Iref* ∈ [0.5, 1.6]A and *Vin* = 10V fixed.

16.5 17.0

16.0

15.5 15.0

*vC(V)*

13.5

0.65

0.50 0.45

iL(A)

0.45

(down) temporal waveform of the current *iL*; (b) phase plane (*vC*, *iL*).

0.50

0.55

0.60

0.65

0.70

0.75

0.55

*iL(A)*

0.60

0.70

14.0 14.5

8.0e−005 1.2e−004

1.2e−004

8.0e−005 1.6e−004

0.0e+000 4.0e−005 2.0e−004

*t(s)*

vC(V) 13.5 14.0 14.5 15.0 15.5 16.0 16.5

(b)

**Figure 8.** Fundamental periodic regime for *Iref* = 0.7A: (a) (up) temporal waveform of the voltage *vC*,

(a)

0.0001

0.0001

0.0001353

Simulation of Piecewise Hybrid Dynamical Systems in Matlab 23

0.0e+000 4.0e−005 1.6e−004 2.0e−004

0.0000353 0.0001353

0.0000353

**Figure 6.** Bifurcation diagram of the Boost in the space (*Iref* , *iL*, *vc*) for *Iref* ∈ [0.5, 1.6]A and *Vin* = 10V fixed.

**Figure 7.** Zoom of figure 5: Border collision for *Iref* = 1.23A.

20 Will-be-set-by-IN-TECH

0.6 0.8 <sup>1</sup> 1.2 1.4 1.6

Iref(A)

**Figure 6.** Bifurcation diagram of the Boost in the space (*Iref* , *iL*, *vc*) for *Iref* ∈ [0.5, 1.6]A and *Vin* = 10V

1.2 1.23 1.4

Iref(A)

0.4 0.6 0.8 1 1.2 1.4 1.6 10 15 20 25 30

iL(A)

0.5 0.6 0.7 0.8 0.9 1 1.1 1.2 1.3 1.4

**Figure 7.** Zoom of figure 5: Border collision for *Iref* = 1.23A.

iL(A)

fixed.

vC(V)

**Figure 8.** Fundamental periodic regime for *Iref* = 0.7A: (a) (up) temporal waveform of the voltage *vC*, (down) temporal waveform of the current *iL*; (b) phase plane (*vC*, *iL*).

t(s)

(a)

vC(V)

(b)

**Figure 10.** Cycle of order 4 for *Iref* = 1.3A: (a) (up) temporal waveform of the voltage *vC*, (down)

12 14 16 18 20 22 24 26

0.0000 0.0001 0.0002 0.0003 0.0004 0.0005 0.0006 0.0007 0.0008

0.0004

0.0006

0.0007114

0.0000 0.0001 0.0002 0.0003 0.0004 0.0005 0.0006 0.0007 0.0008

0.0003114

0.0000894 0.0004894

0.0001065 0.0003114 0.0005065 0.0007114

0.0004 0.0006

Simulation of Piecewise Hybrid Dynamical Systems in Matlab 25

0.0005 0.0005065

0.0004894

4T

0.0002

0.0001065 0.0001

0.0000894

0.0002

0.0001 0.0005

vC(V)

iL(A)

1.2 1.3 1.4

0.5

iL(A)

0.5

temporal waveform of the current *iL*; (b) phase plane (*vC*, *iL*).

0.6

0.7

0.8

0.9

1.0

1.1

1.2

1.3

1.4

0.6 0.7 0.8 0.9 1.0 1.1

**Figure 9.** Cycle of order 2 for *Iref* = 1A:(a) (up) temporal waveform of the voltage *vC*, (down) temporal waveform of the current *iL*; (b) phase plane (*vC*, *iL*).

22 Will-be-set-by-IN-TECH

0.0002

0.0001 0.0002

t(s)

vC(V)

(b)

**Figure 9.** Cycle of order 2 for *Iref* = 1A:(a) (up) temporal waveform of the voltage *vC*, (down) temporal

13 14 15 16 17 18 19 20 21

(a)

0.0000 0.0001 0.0002 0.0003 0.0004

0.0002

2T

0.0000 0.0003

0.0001096

0.0001

0.0000805

0.0000805 0.0001096 0.0001

0.0004

0.0003096

0.0003

0.0002804

0.0002804 0.0003096 0.0003

18 19

20

21

17

vC(V)

iL(A)

16

14 15

13

1.0

1.1

0.7 0.6 0.5

0.4

iL(A)

0.4

waveform of the current *iL*; (b) phase plane (*vC*, *iL*).

0.5

0.6

0.7

0.8

0.9

1.0

1.1

0.8 0.9

**Figure 10.** Cycle of order 4 for *Iref* = 1.3A: (a) (up) temporal waveform of the voltage *vC*, (down) temporal waveform of the current *iL*; (b) phase plane (*vC*, *iL*).

24 Will-be-set-by-IN-TECH 26 MATLAB – A Fundamental Tool for Scienti c Computing and Engineering Applications – Volume 3 Simulation of Piecewise Hybrid Dynamical Systems in Matlab <sup>25</sup>

**4. Conclusion**

current mode.

**Acknowledgments**

**Author details** Fatima El Guezar

Hassane Bouzahir

**5. References**

633-647.

E–14.

pp. 359-371.

Catalunya, Spain.

Appliquées de Toulouse, France.

*Electronic Systems*, Tokushima, Japan, pp. 62–65.

In this chapter, we have showed an accurate and fast method to determine events' occurrence for planar piece-wise affine hybrid systems. As a result, we have implemented our algorithm

Simulation of Piecewise Hybrid Dynamical Systems in Matlab 27

This toolbox has also been completed by analysis tools such as displaying the bifurcation and parametric diagrams. The algorithm takes the advantage of the analytical form that appears in the planar case. Our approach can not be extended to a higher dimension. DC-DC converters like Boost converter are known to be simple switched circuits but very rich in nonlinear dynamics. As application, we have chosen the example of Boost converter controlled in

The authors would like to thank Pascal Acco and Danièle Fournier–Prunaret for crucial

discussions on the original version of our work on this subject.

*ESSI & ERMAGIM, Ibn Zohr University, EST, PO Box 32/S, Agadir, Morocco*

*ESSI & ERMAGIM, Ibn Zohr University, EST, PO Box 32/S, Agadir, Morocco*

*Faculty of Engineering, AlHosn University, PO Box 38772, Abu Dhabi, United Arab Emirates*

[1] Acco, P. (December 2003). Etude de la boucle à verrouillage de phase par impulsions de charge: Prise en compte des aspects hybrides. *Ph D thesis*, Institut National des Sciences

[2] Banerjee, S. (2000). Bifurcations in two-dimensional piecewise smooth maps - theory and applications in switching circuits, *IEEE Trans. on Circuits and Systems-I*, Vol.47, pp.

[3] Bouzahir H.; El Guezar, F. & Ueta, T. (2007). On Scicos simulation of a hybrid dynamical system. *Proceedings of the 15th IEEE International Workshop on Nonlinear Dynamics of*

[4] Brockett R. W. & Wood J. R. (1984). Understanding power converter chaotic behavior mechanisms in protective and abnormal modes. *Proceedings of POWERCON 11*, pp.

[5] Deane, J. H. B. & Hamill D. C. (1990). Instability, subharmonics, and chaos in power electronic systems. *IEEE Trans. Power Electronics*, Volume 5, pp. 260–268, 1990. [6] El Aroudi, A.; Benadero, L.; Toribio, E. & Machiche, S. (2000). Quasiperiodicty and chaos in the DC-DC buck-boost converter. *International Journal of Bifurcation and Chaos*, Vol. 10,

[7] El Aroudi, A. (February 2000). Study of nonlinear phenomena and quasiperiodicity route to chaos in PWM DC/DC converters. *Ph D thesis*, Universitat Politécnica de

in Matlab toolbox version (free downloadable on http://felguezar.000space.com/).

**Figure 11.** Chaotic regime for *Iref* = 1.5A: (a) (up) temporal waveform of the voltage *vC*, (down) temporal waveform of the current *iL*; (b) phase plane (*vC*, *iL*).

## **4. Conclusion**

24 Will-be-set-by-IN-TECH

t(s)

vC(V)

30

10 15 20 25

**Figure 11.** Chaotic regime for *Iref* = 1.5A: (a) (up) temporal waveform of the voltage *vC*, (down)

(b)

(a)

0.0000 0.0005 0.0010 0.0015 0.0020 0.0025 0.0030 0.0035 0.0040

0.0000 0.0005 0.0010 0.0015 0.0020 0.0025 0.0030 0.0035 0.0040

iL(A)

vC(V)

10

1.6

0.6 0.7 0.8 0.9 1.0 1.1 1.2 1.3 1.4 1.5

0.5

iL(A)

0.5

temporal waveform of the current *iL*; (b) phase plane (*vC*, *iL*).

0.6

0.7

0.8

0.9 1.0 1.1 1.2 1.3 1.4 1.5 1.6

15

20

25

30

In this chapter, we have showed an accurate and fast method to determine events' occurrence for planar piece-wise affine hybrid systems. As a result, we have implemented our algorithm in Matlab toolbox version (free downloadable on http://felguezar.000space.com/).

This toolbox has also been completed by analysis tools such as displaying the bifurcation and parametric diagrams. The algorithm takes the advantage of the analytical form that appears in the planar case. Our approach can not be extended to a higher dimension. DC-DC converters like Boost converter are known to be simple switched circuits but very rich in nonlinear dynamics. As application, we have chosen the example of Boost converter controlled in current mode.

## **Acknowledgments**

The authors would like to thank Pascal Acco and Danièle Fournier–Prunaret for crucial discussions on the original version of our work on this subject.

## **Author details**

Fatima El Guezar *ESSI & ERMAGIM, Ibn Zohr University, EST, PO Box 32/S, Agadir, Morocco*

Hassane Bouzahir

*ESSI & ERMAGIM, Ibn Zohr University, EST, PO Box 32/S, Agadir, Morocco Faculty of Engineering, AlHosn University, PO Box 38772, Abu Dhabi, United Arab Emirates*

## **5. References**

	- [8] El Aroudi, A. & Leyva, R. (2001). Quasi-periodic route to chaos in a PWM voltage-controlled DC-DC boost converter. *IEEE Trans. on Circuits and Systems*, Vol. 48, No. 8, pp. 967-978.
	- [9] El Aroudi, A.; Debbat, M.; Giral, R.; Olivar, G.; Benadero, L.& Toribio, E. (2005). Bifurcations in DC-DC switching converters: review of methods and applications. *International Journal of Bifurcation and Chaos*, Vol. 5, pp. 1549–1578.
	- [10] El Guezar, F. & Bouzahir, H. (2008). Chaotic behavior in a switched dynamical system. *Modelling and Simulation in Engineering*, Vol. 2008, Article ID 798395, 6 pages.
	- [11] El Guezar, F.; Acco, P.; Bouzahir, H. & Fournier-Prunaret, D. (2008). Accurate and Fast Event Detection Occurrence in Planar Piecewise Affine Hybrid Systems. *Proceedings of the International Symposium NOLTA (NOn Linear Theory and its Applications)*, September 7-10, Budapest-Hungary, pp. 341–344.
	- [12] El Guezar, F. (December 2009). Modélisation et simulation des systèmes dynamiques hybrides affines par morceaux. Exemples en électronique de puissance. *Ph D thesis*, Institut National des Sciences Appliquées de Toulouse, France.
	- [13] El Guezar, F.; Bouzahir, H. & Fournier-Prunaret, D. (2011). Event Detection Occurrence For Planar Piecewise Affine Hybrid Systems. *Nonlinear Analysis: Hybrid Systems*, Vol. 5, pp. 626–638.
	- [14] Girard, A. (2002). Detection of event occurrence in piece-wise linear hybrid systems. *Proceedings of the 4th International Conference on Recent Advances in Soft Computing*, Nottingham, United Kingdom, December, pp. 19–25.
	- [15] Girard, A. (September 2004). Analyse algorithmique des systèmes hybrides. *Ph D thesis*, Institut National Polytechnique de Grenoble, France.
	- [16] Hamill , D. C. & Jeffries, D. J. (1988). Subharmonics and chaos in a controlled switched-mode power converter. *IEEE Trans. Circuits Systs. I*, Vol. 35, No. 8, pp. 1059–1060.
	- [17] Hedayat, C. D.; Hachem, A.; Leduc, Y. & Benbassat, G. (March–April 1997). High-level modeling applied to the second-order charge-pump PLL circuit. *Technical report, Texas Instrument Technical Journal*. Vol. 14, No. 2.
	- [18] Mira, C. (1987). Chaotic dynamics. *World scientific Publishing*.
	- [19] Tse, C.K. (1994). Chaos from a Buck switching regulator operating in discontinuous mode. *IEEE Transactions on International Journal of Circuit Theory and Application*. Vol. 22, No. 7, pp. 262–278.
	- [20] Tse, C.K. (1994). Flip bifurcation and chaos in three–state boost switching regulators. *IEEE Transactions on Circuits and Systems I: Theory and Applications*, Vol. 41, No. 1, pp. 16–23.
	- [21] Tse, C.K. (2003). *Complex behavior of switching Power converters*, CRC Press.
	- [22] Van Paemel, M. (July 1994). Analysis of a charge pump PLL: a new model. *IEEE Transactions on Communications*, Vol. 42, No. 7, pp. 2490-2498.
	- [23] Yuan, G. H.; Banerjee, S.; Ott, E. & Yorke, J. A. (1998). Border-collision bifurcations in the buck converter, *IEEE Trans. on Circuits and Systems-I*, Vol. 45, pp. 707–715.

© 2012 Belavý et al., licensee InTech. This is an open access chapter distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

© 2012 Belavý et al., licensee InTech. This is a paper distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

**Robust Control of Distributed Parameter** 

**Systems with Demonstration in Casting** 

Cyril Belavý, Gabriel Hulkó and Karol Ondrejkovič

Additional information is available at the end of the chapter

**Software Support** 

http://dx.doi.org/10.5772/46460

**1. Introduction** 

**Technology and MATLAB/Simulink/DPS Blockset** 

Most of the dynamical systems analysed in engineering practice have the dynamics, which depends on both position and time. Such systems are classified as distributed parameter systems (DPS). The time-space coupled nature of the DPS is usually mathematically described by partial differential equations (PDE) as infinite-dimensional systems. However, from point of view of implementation of DPS control in technological practice, where a finite number of sensors and actuators for practical sensing and control is at disposal, such infinite-dimensional systems need to be approximated by finite-dimensional systems. There

In the first mathematical foundations of DPS control, analytical solutions of the underlying PDE have been used (Butkovskij, 1965; Lions, 1971; Wang, 1964). That is the decomposition of dynamics into time and space components based on the eigenfunctions of the PDE. Continuous and approximation theories aimed to control of parabolic systems presents monograph (Lasiecka & Triggiani, 2000). Methodical approach from the view of time-space separation with model reduction is presented in (Li & Qi, 2010). Variety of transfer functions for systems described by PDE are illustrated by means of several examples in (Curtain & Morris, 2009). Well-known reduction methods based on finite difference method (FDM), or finite element method (FEM), spectral method require an accurate nominal PDE model and

are many dimension reduction methods, which can be used to solve this problem.

usually lead to a high-order model, which requires unpractical high-order controller.

An engineering approach for the control of DPS is being developed since the eighties of the last century (Hulkó et al., 1981, 1987, 1998, 2009a, 2009b). In the field of lumped parameters system (LPS) control, where the state/output quantities *x(t)/y(t)* – parameters are given as
