Sample-and-hold control

Consider the following sample-and-hold architecture:

Image

where the sampling time is not constant. A new sample is obtained when \left|m-z\right|\geqslant\delta, \delta>0. The memory m remains constant until a new sample is obtained.

In order to model the closed loop model as a hybrid system, we must define C,f,D,g

State Vector:

x=\left[\begin{array}{c}  z\\m\\u\end{array}\right]

f(x)=\left[\begin{array}{c}  Az\\ 0\\ 0\\\end{array}\right]

C=\left\{ \mathbb{R}^{n_{p}+n_{c}+n_{p}}\,:\left|m-z\right|<\delta\right\}

g(x)=\left[\begin{array}{c}  z\\-Kz\\-Kz\end{array}\right]

D=\left\{ \mathbb{R}^{n_{p}+n_{c}+n_{p}}\,:\left|m-z\right|\geqslant\delta\right\}

The simulation results are shown in the following images:

Image

Image

Image

Simulation code

%file name: run.m
clear all;
close all;
%%
% planta con S&H, se guarda una muestra de la planta en "m".
% No se tomara otra muestra hasta que la distancia entre la salida de la
% planta y m es mayor que DELTA

%planta:
A = -.5; %B=0;
%realimentación
K = 1; DELTA = 3;

%% Define initial conditions:
x0 = [10; 0; -10];

%% Simulation horizon
TSPAN = [0 10];
JSPAN = [0 15];
rule = 1;
options = odeset('RelTol', 1e-6, 'MaxStep', .1);

[t j x] = HyEQsolver(@f,@g,@C,@D,x0,TSPAN,JSPAN,rule,options);

%%
plotHybridArc(t,j,x(:,3))
title('Input')
xlabel('j')
ylabel('t')
zlabel('entrada')

plotHybridArc(t,j,abs(x(:,2)-x(:,1)))
title('Delta')
xlabel('j')
ylabel('t')
zlabel('delta ')

plotHybridArc(t,j,x(:,1))
title('Output')
xlabel('j')
ylabel('t')
zlabel('z')

flow map

%file name: f.m
 function [ xdot ] = f( x )
 %F Summary of this function goes here
 % Detailed explanation goes here

A = -0.5; %B=0;
 K = 1; DELTA = 3;

z = x(1);
 zdot = A*z;

xdot = [zdot;0;0];
 end

jump map

function [ xPlus] = g(x)
%UNTITLED Summary of this function goes here
% Detailed explanation goes here

A = -.5; %B=0;
K = 1; DELTA = 3;

z = x(1);
zPlus = z;
mPlus = z;
uPlus = -K*z;

xPlus = [zPlus; mPlus; uPlus];
end

flow condition

function [ inC ] = C( x)
%INC Summary of this function goes here
% Detailed explanation goes here

A = -0.5; %B=0;
K = 1; DELTA = 3;
%%%

z = x(1);
m = x(2);
if abs(z-m)<DELTA
inC = true;
else
inC = false;
end

end

jump condition

function [ inD ] = D( x)
%INC Summary of this function goes here
% Detailed explanation goes here

A = -0.5; %B=0;
K = 1; DELTA = 3;
%%%

z = x(1);
m = x(2);
if abs(z-m)>=DELTA
inD = true;
else
inD = false;
end

end

Advertisements

Tire Pressure Control

Model

Let’s consider a pneumatic system, which consists of a flat tire and an air compressor to inflate it. The system can be modeled as follows:

Modelo
where:

R_c : is the compressor’s pneumatic resistance.

C_c : models the tire’s capacitance

R_p : models the losses due to the puncture.

To develop a hybrid model I chose the following variables:

p : stands for the manometric pressure in the tire.

q : stands for the state of the compressor (‘on’ or ‘off’).

The compressor turns on if the pressure falls below p_{min} and it turns off when the pressure inside the tire reaches p_{max} . Between this two limits the pressure flows according to the following equation:

\dot{p}=\left(-\frac{R_{c}+R_{p}}{R_{c}R_{p}C_{c}}p+\frac{P_{comp}}{R_{c}C_{c}}\right)q+\left(-\frac{C_{c}}{R_{p}}p\right)\left(1-q\right)

This equation was obtained under the following assumptions:

• The variation of the tire’s volume is negligible, then tire’s capacitance is constant.

• The compressor is a constant pressure source of value P_{comp}

• The pressure never exceeds the limits p_{min}  or p_{max}

The hybrid model is then defined by:

\mathcal{H}=(C,f,D,g)

x=\left[\begin{array}{c}  p\\  q  \end{array}\right]

C=\{x\epsilon\mathbb{R}^{2}/p_{min}+\varepsilon\leq p\leq p_{max}-\varepsilon\}

D=\{x\epsilon\mathbb{R}^{2}/\left|p-p_{min}\right|<\varepsilon\,\vee\,\left|p-p_{max}\right|<\varepsilon\}

f(x)=\left[\begin{array}{c}  \left(-\frac{R_{c}+R_{p}}{R_{c}R_{p}C_{c}}p+\frac{P_{source}}{R_{c}C_{c}}\right)q+\left(-\frac{C_{c}}{R_{p}}p\right)\left(1-q\right)\\  0  \end{array}\right]

g(x)=\left\{ \begin{array}{c}  \left[\begin{array}{c}  11-\varepsilon\\  0  \end{array}\right]\,,\, if\,\left|p-p_{max}\right|<\varepsilon\\  \left[\begin{array}{c}  8+\varepsilon\\  1  \end{array}\right]\,,\, if\,\left|p-p_{min}\right|<\varepsilon  \end{array}\right.

Simulation

The Matlab code written to simulate the system using HyEQ Toolbox is presented below. The values of the parameters were chosen in order to obtain results that could be quickly visualized,  but they must be changed for more appropiate values.

C.m

function inC = C(x)
%check if in C

p = x(1);
eps = 0.12;
pmin = 8;
pmax = 11;

if (p>=pmin+eps && p<=pmax-eps) || p>=pmax+eps
 inC = 1;
else
 inC = 0;
end

end

D.m

function inD = D(x)
%check if in D

p = x(1);
eps = 0.12;
eps = 0.12;
pmin = 8;
pmax = 11;
if (p<=pmin+eps && p>=pmin-eps ) || (p<=pmax+eps && p>=pmax-eps )
 inD = 1;
else
 inD = 0;
end
end

f.m

function xdot = f(x)
%flow map

p = x(1);
q = x(2);
if q == 1
 %pdot=-(Rc+Rp)/(Rc*Rp*Cc)*p+psource/(Cc*Rc)
 pdot = -0.2*p+0.3*7.5;
else
 %pdot=-Cc/Rp*p
 pdot = -0.1*p;
end
qdot = 0;

xdot = [pdot qdot]';
end

g.m

function xplus = g(x)
% jump map

p = x(1);
q = x(2);

pmin = 8;
pmax = 11;
eps = 0.12;

if p<=pmax+eps && p>=pmax-eps
 qplus = 0;
 pplus = pmax-eps;
elseif p<=pmin+eps && p>=pmin-eps
 qplus = 1;
 pplus = pmin+eps;
end

xplus = [pplus qplus]';

end

Results

sim1 sim2

Hybrid model for a switching system

Consider a system with a dc motor and the load, coupled by flexible shaft as shown in the next diagram

motorresortemasa

where:
\varphi \; \dot{\varphi} : Position and speed of the load respectively.

\theta \; \dot{\theta} : Position and speed of the motor shaft respectively.

B_{L}\;B_{m}: Dynamic friction coefficient of the load and the motor respectively.

J_{L}\;J_{m}: Inertia moment of the load and the motor respectively.

K_{s}: Elastic torsion coefficient of the spring.

The states of system are the position and speed of the motor and the load respectively, z=[\varphi \; \dot{\varphi} \; \theta \; \dot{\theta}] and the input is the voltage applied. With this and using Newtons Law and the electrical equations for the motor, the system in state space is:

\left\{ \begin{array}{ccc}  \dot{z} & = & Az+Bv\\  \varphi & = & Cz  \end{array}  \right.

\mbox{A}= \left( \begin{array}{cccc}  0 & 1 & 0 & 0\\  -\frac{K_{s}}{J_{L}} & -\frac{B_{L}}{J_{L}} & \frac{K_{s}}{J_{L}} & 0\\  0 & 0 & 0 & 1\\  \frac{K_{s}}{J_{M}} & 0 &-\frac{K_{s}}{J_{M}} & -\frac{K^{2}+R_{A}\;B_{M}}{J_{M}\;R_{A}}  \end{array}  \right) \mbox{B}=\left( \begin{array}{c}  0\\  0\\  0\\  \frac{K}{J_{M}\;R_{A}}\\  \end{array}  \right) \mbox{C}=\left( \begin{array}{lcccr}  1 & 0 & 0 & 0\\  \end{array}  \right)

To control the position in a point of reference a sliding mode control with next structure is developed

esquema_lineal_md

Because of the implemented control the signal input to plant is not continuous, and it can be modelled as a switching signal with hysteresis, that jumps between two differences values depending of signal e .
The close loop system can be modelled as a hybrid system, if an extra variable \tau which count the time between jumps is considered and is reseted when jump occurs. the variable that jumps is called q .

x = \left[ \begin{array}{c}  z\\  \tau\\  q  \end{array}  \right]

\dot{x} = \left[ \begin{array}{c}  f_{q}\\  \dot{\tau}\\  \dot{q}  \end{array}  \right] = \left[ \begin{array}{c}  Ax+B(q\,\mbox{v})\\  1\\  0  \end{array}  \right] \qquad C=\{ z \in\mathbb{R}^{4} \; , \tau \in \mathbb{R}: \tau \ge 0 \; , q \in \mathbb{Q} \}

x^{+} = \left[ \begin{array}{c}  z^{+}\\  \tau^{+}\\  q^{+}  \end{array}  \right] = \left[ \begin{array}{c}  z\\  0\\  -q  \end{array}  \right]

D = \{(q,u) \in \{-1;1\} \times \mathbb{R}: q=-1 , u \ge \overline{h} \} \; \bigcup \; \{ (q,u) \in \{-1;1\} \times \mathbb{R}: q=1 , u \leq \underline{h} \}

Where \overline{h} ,\; \underline{h} are the levels of hystereses and v is the positive bounded of the switching signal.

solutionarc

% run.m to simulate timer with reset

clear all

global u
global r

%Reference Position [rad]
r=5;

%Initial Value for input for de switching signal
u=0;

%Define inicial conditions of the states

x0=[0;0;0;0;0;1];

%simulation horizon
TSPAN = [0 4];
JSPAN = [0 15];

%ruls for jump
%rule  = 1 -> priority for jumps
%rule  = 2 -> priority for flows
%rule  = 3 -> ramdom priority

rule=1;

options = odeset('RelTol',1e-6,'MaxStep',.1);

%run simulation
[t j x] = HyEQsolver(@f,@g,@C,@D,x0,TSPAN,JSPAN,rule,options);

figure(1)
clf
subplot(2,1,1),plotflows(t,j,x(:,1))
grid on

ylabel('x1')
subplot(2,1,2),plotjumps(t,j,x(:,1))
grid on
ylabel('x1')

% plot hybrid arc
plotHybridArc(t,j,x)
xlabel('j')
ylabel('t')
zlabel('x')
function inD = D(x)
% Jump set
% Check if in D

global u

% Hysterisis level
hbar=1;

q=x(6);

if ((q == -1) && (u >= hbar))
    inD = 1;
elseif ((q==1) && (u <= -hbar))
    inD = 1;
else
    inD=0;
end

function inC = C(x)
%flow set
%checck if in C

  inC=1;

end

function xdot = f(x)
%flow map

global u
global r

z1=x(1);
z2=x(2);
z3=x(3);
z4=x(4);

tau=x(5);
q=x(6);

z=[z1;z2;z3;z4];

% Linear System
A=[0 1 0 0 ; -90.8 -0.7 90.8 0 ; 0 0 0 1 ; 55,41 0 -55.41 - 105.6];
B=[0 ; 0 ; 0 ; 220.7];

K=[2.2895 0.2714 0.8352 0.0045];

u = r-(K*z);

%Flow
xdot = [(A*z)+(B*(q*24)) ; 1 ; 0];

end

Hybrid Brain Waves

Visual Spatial Covert Attention is a cognitive effect that is produced when a person is looking into whatever it has in front of him, and some event appears on the visual field.  If the person decides to pay attention to the event without gazing at it,  a shift is produced in the  alpha band of the EEG signals (10 Hz) towards those electrodes located on the same side of the event location.

Visual Spatial Covert Attention effect

Visual Spatial Covert Attention effect

 

 

Events (q) \  are triggered one per second according to the timer (\tau) \  on each side of the screen, one after the other, being -1 for left, and 1 for right.  For each event, the brain signals will shift towards each side (u) \ and if a threshold is reached, the output of the system (y) \ will change to the side where the brain signals drifted.

 

 

 

EEG Channels

EEG Channels

 

 

Each signal from each channel is band-passed to 8-12 hz and the power spectral density of all of them will be  vectorized and summed up (according to their relative position, as suggested by the the figure)

We are considering that the person always pays COVERT attention to the event, and the signal will drift to the same side where each event was triggered.

 

 

 

The model is: \mathcal{H}=(C,f,D,g) \

x =\left[  \begin{array}{c}  u\\  y\\  \tau\\    q\end{array}    \right]

f(x) =\left[    \begin{array}{c}    p(x) \\    0\\    1\\    0\\    \end{array}  \right]

For this model, we will consider the brainwaves to be linear drifts to each direction.

p(x) =\left[  \begin{array}{ll}  1 & q = 1 \\    -1 & q = -1 \\    \end{array}  \right]

What about the jump-map  ?

g(x) =\left[  \begin{array}{c}  0 \\  s(x) \\  0\\  -q\\  \end{array}  \right]

s(x) =\left[  \begin{array}{ll}  1 & u \ge h_{lower}\\    -1 & u \leq h_{upper} \\    \end{array}  \right]

The sets being:

C:=\{ x \in \mathbf{R} \times \mathbf{R} \times [0,1] \times {0,1} \}

D:=\{ x \in \mathbf{R} \times \mathbf{R} \times [0,1] \times {0,1} : \tau = 1 \} \cup \{ x \in \mathbf{R} \times \mathbf{R} \times [0,1] \times {0,1}: y=-1, u \ge h_{upper} \} \cup \{ x \in \mathbf{R} \times \mathbf{R} \times [0,1] \times {0,1}: y=1, u \le h_{lower} \}

Ok, let’s see the results on the simulation:

System's Dynamics Simulation

System’s Dynamics Simulation

It can be seen how the output of the system, adapts accordingly to the direction where the brainwaves drifted (left, right) based on the direction of the previous stimulation event.   The timer is only triggered the first time because the brainwaves are moving towards the stimulus’ direction (left) and the output is already there.

The output of the system is driven by the drift of the brainwaves.

The output of the system is driven by the drift of the brainwaves.

 

 

If we plot u against y  we see a sheared hysteresis characteristic plot.  The brainwaves affects the output of the system after the jump has occurred, so both signals end up being shifted by a lag.

 

In red, jumps are plotted from D   regions to their new positions after the jump.

 

 

 

Simulation Code:

function inC = C(x)
% check if x is in C
tau = x(3);
if  (tau >= 0) && (tau <= 1)
inC = 1;
else
inC = 0;
end
end
function inD = D(x)
% check if x is in D
global h_upper;
global h_lower;
tau = x(3);
u = x(1);
y = x(2);
if  (tau >= 1) || (y == -1 && u>=h_upper) || (y==1 && u<=h_lower)
      inD = 1;
else
inD = 0;
end

end
function xdot = f(x)
% flow map

q=x(4);
y=x(2);

if (q==1)
 udot = 1;
else
 udot = -1;
end

xdot = [ udot; 0; 1; 0];

end

 

function xplus = g(x)
% jump map

global h_upper;
global h_lower;
u = x(1);
q = x(4);

% on jumps we are always over or bellow the bound.
if (u>=h_upper)
 ydot = 1;
elseif (u<=h_lower)
 ydot = -1;
else
 ydot = x(2);
end
xplus = [0; ydot; 0; -q];

end

Finally, the run.m code


function hybridsystem()

global h_upper;

global h_lower;

h_upper = 0.5;

h_lower = -0.5;

% define initial conditions

x0 = [0; -1; 0; -1]

% simulation horizon

TSPAN = [0 10];

JSPAN = [0 20];

% Rule for OVERLAPS, what to choose when C and D collide?

% rule = 1 -> priority for jumps

% rule = 2 -> priority for flows,

% rule = 3 -> random priority

rule=1;

options = odeset('RelTol', 1e-6,'MaxStep', .1);

%Do your magic.....

[t j x] = HyEQsolver( @f,@g,@C,@D,x0, TSPAN, JSPAN, rule, options);

% t continuous time

% j  discrete time

% x state variable for each (j,t)

figure(1);

clf

subplot(4,1,1);plotflows(t,j,x(:,1));

grid on;

ylabel('u');

&nbsp;

subplot(4,1,2);plotflows(t,j,x(:,2));

grid on;

axis([0 10 -1.5 1.5]);

ylabel('y');

&nbsp;

subplot(4,1,3);plotflows(t,j,x(:,3));

grid on;

ylabel('tau');

&nbsp;

subplot(4,1,4);plotflows(t,j,x(:,4));

grid on;

axis([0 10 -1.5 1.5]);

ylabel('q');

&nbsp;

figure(2);

plot(x(:,1),x(:,2),'-');

grid on;

ylabel('y: output');

xlabel('u: brain waves');

axis([-2 2 -2 2]);

Done!.  Thanks for reading !

Temperature Control of an Alkaline Electrolyser

The chosen Hybrid system is a model aproximation to a temperature control used at ITBA’s Alkaline Electrolyser. Electrolysis is an exothermic process which produces heating of KOH solution. Intending to analyze process efficience as a function of temperature, it is considered a group of temperatures to make tests. With this results, relationships between power consumed, production of gases and their purity will be concluded.
To maintain system temperature T_1 , it is considered natural convection \dot{Q}_{NC} with environment temperature T_0 of the entire equipment, and dissipated power by a cooling system composed of a countercurrent flow of water and a radiator which also interacts with the environment, \dot{Q}_{Cool} . For this approximation it will be considered without the intermediate step of cooling water: \dot{Q}_{Cool} = K_{Cool}.(T_0 - T_1) .
This system can be considered Hybrid because of the Cooling System Switching. Through state of q, it is shown the activation of \dot{Q}_{Cool} .

Then, the proposed model is:
\mathcal{H}=(C,f,D,g) \\ \\  x=\left[ \begin{array}{c} x \\ q \end{array} \right] \\ \\  C=\left\{x\in\mathbb{R}^2/(T > T_{min}\; \& \;q=0)\;|\; (T < T_{MAX}\; \& \;q=1) \right \} \\ \\  D=\left\{x\in\mathbb{R}^2/(T\geqslant T_{MAX}\; \& \;q=0)\;|\; (T\leq T_{min}\; \& \;q=1) \right \} \\ \\  f(x)=\begin{bmatrix} 0 \\ \frac{\dot{Q}+[K_{NC}+q.K_{Cool}.(T_0-T_1)]}{C} \end{bmatrix} \\ \\ \\  g(x)=\begin{bmatrix} 1-q \\ T \end{bmatrix}

Used parameters are:
T_{Min} = 40 \;^{o}\textrm{C} \\  T_{Max} = 45 \;^{o}\textrm{C} \\  T_{Env} = 20 \;^{o}\textrm{C} \\  \dot{Q} = 1000 \;kW \\  K_{NC} = 5 \;kW/K \\  K_{Cool} = 50 \;kW/K

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Function run
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function run

Parameters

% initial conditions
q_0 = 0;
Temp_0 = 20;
x0 = [q_0;Temp_0];

% simulation horizon
TSPAN=[0 1000];
JSPAN = [0 20];

% rule for jumps
% rule = 1 -> priority for jumps
% rule = 2 -> priority for flows
rule = 1;

options = odeset('RelTol',1e-6,'MaxStep',.1);

maxStepCoefficient = .1;  % set the maximum step length. At each run of the
                   % integrator the option 'MaxStep' is set to
                   % (time length of last integration)*maxStepCoefficient.
                   %  Default value = 0.1

% simulate
[t x j] = HyEQsolver( @f,@g,@C,@D,x0,TSPAN,JSPAN,rule,options,maxStepCoefficient);

% plot solution
figure(1) % position
clf
subplot(2,1,1),plotflows(t,j,x(:,1))
grid on
ylabel('q')

subplot(2,1,2),plotjumps(t,j,x(:,1))
grid on
ylabel('q')

figure(2) % velocity
clf
subplot(2,1,1),plotflows(t,j,x(:,2))
grid on
ylabel('Temp')

subplot(2,1,2),plotjumps(t,j,x(:,2))
grid on
ylabel('Temp')

% plot hybrid arc
plotHybridArc(t,j,x(:,2))
xlabel('j')
ylabel('t')
zlabel('Temp')

% figure(3)
% plot(t,y(:,1));

Parameters (Parameters.m)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Definition of Parameters for Model of Cooling Hysteresis
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

global TempMin TempMax Q Knatconv Kcool TempEnv

TempMin = 50;   % °C
TempMax = 55;   % °C

Q = 1000;       % kW
Knatconv = 5;   % kW/K
Kcool = 50;   % kW/K
TempEnv = 20;   % °C

Flow map (f.m)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Declaration of Flow Map F
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function xdot = f(x)

% Parameters
global Q        % kW
global Knatconv	% kW/K
global Kcool	% kW/K
global TempEnv	% °C

% State variables
q = x(1);
Temp = x(2);

% Differential equations
xdot = [0 ; (Q + (Knatconv + q * Kcool) * (TempEnv - Temp)) / 1000];
end

Flow set (C.m)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Declaration of Flow Set C
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function inC = C(x)

% Out:
%   0 if value is outside C
%   1 if value is inside C

global TempMin
global TempMax

q = x(1);
Temp = x(2);

if (((Temp < TempMax) && (q == 0)) || ((Temp > TempMin) && (q == 1)))
    inC = 1;
else
    inC = 0;
end
end

Jump map (g.m)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Declaration of Jump Map G
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function xplus = g(x)

% State variables
q = x(1);
Temp = x(2);

xplus = [1-q ; Temp];
end

Jump set (D.m)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Declaration of Jump Set D
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function inD = D(x)

% Out:
%   0 if value is outside D
%   1 if value is inside D

global TempMin
global TempMax

q = x(1);
Temp = x(2);

if (((Temp <= TempMin) && (q == 1)) || ((Temp >= TempMax) && (q == 0)))
    inD = 1;
else
    inD = 0;
end
end

Figures 1 and 2 represents temperature T_1 and Hybrid arc \phi , as a function of continuous and discrete time, t and j .

Fig 1. Temp VS t-j

Fig 2. Temp VS arc

Modelo Híbrido del SQUID

Graduate Course on Robust Hybrid Control Systems – ITBA

Homework 1

Se propone y simula un modelo híbrido para un sensor SQUID (Superconducting Quantum Interference Device) con realimentación de flujo magnético.

Este sensor es un dispositivo no lineal que utiliza un lazo de realimentación para incrementar su rango dinámico. Suele ser modelado como un filtro pasa bajos con ganancia y frecuencia de corte conocida. Señales del alto slew rate pueden producir corrimientos en el punto de trabajo y generar discontinuidades en la tensión de salida.

Como modelo híbrido se propone un sistema conmutado que incluye la señal de entrada al sensor, la dinámica del mismo y se asume una distribución de probabilidad del nuevo punto de trabajo al realizarse un salto.

El vector de estados se define como una variable continua Z \in [-10,10] que representa la salida del SQUID, un contador \tau \in [0,\tau_D], un indicador del punto de trabajo q\in \mathcal{Z}, la entrada del sistema \phi_1 \in \Re y su derivada \phi_2 \in \Re. En este trabajo se estudian señales triangulares con valor pico \phi_{pico} y período (2\cdot \tau_D).

Los saltos se realizan de manera periódica en un cuarto del período de la señal de entrada. En cada salto, el contador se reinicia en cero, la variable continua se mantienen igual, se modifica aleatoriamente el punto de trabajo y la señal de entrada se mantiene constante pero se invierte la polaridad de la pendiente. El conjunto y la función de salto quedan definidos como:

eq2Entre los instantes de saltos el sistema fluirá continuamente, la variable Z será modificada según la dinámica del sistema, el punto de trabajo y la señal de entrada. El contador se incrementará en cada instante de simulación mientras que q se mantendrá constante. La señal de entrada modificará su valor linealmente según su pendiente. El conjunto y la función de flujo resultan:

eq3Simulación

A continuación se muestra los resultados de la simulación del sistema propuesto. Para la ganancia en DC y el polo del sensor se tomo S =0.75 V y P=50 kHz. El período  y el valor pico de la señal se tomaron igual a 2 ms y a 5.

3dFigura 1: trayectoria de la variable Z con respecto al tiempo híbrido.

fig2Figura 2: variables Z, \tau y \phi_1 con respecto al tiempo continuo.

fig3Figura 3: ampliación del instante del primer salto en la figura 2 con respecto al tiempo. La dinámica del sensor es visible.

%Flow map
function [xdot]=f(x)
z = x(1);
tau = x(2);
q = x(3);
fi1 = x(4);
fi2 = x(5);

zdot = ( (fi1-q)*0.75 -z)*(2*3.14*50e3);
taudot = 1;
qdot = 0;
fi1dot = fi2;
fi2dot = 0;

xdot = [zdot; taudot; qdot; fi1dot; fi2dot]; 

 


% jump map
function [xplus]=g(x)
z = x(1);
tau = x(2);
q = x(3);
fi1 = x(4);
fi2 = x(5);

qplus=q;
while qplus==q,
 qplus=ceil(random('norm',0,4));
end
zplus=z;
tauplus=0;
fi1plus = fi1;
fi2plus = -fi2;

xplus=[zplus; tauplus; qplus; fi1plus; fi2plus]; 

 


%check if in D
function ind=D(x)
z = x(1);
tau = x(2);
q = x(3);
fi1 = x(4);
fi2 = x(5);
if (tau>=0.001),
 ind=1;
else
 ind=0;
end 

 


%check if in C
function inC=C(x)
z = x(1);
tau = x(2);
q = x(3);
fi1 = x(4);
fi2 = x(5);

if(tau>=0)
 inC=1;
else
 inC=0;
end 

 


%simulate the SQUID
clear all;
close all;

%define inicial condition
T = 0.001; % Período igual a 1 segundo.
x0=[-5*0.75;0;0;-5;10/T]; % Condiciones inicial

%simulation horizon
TSPAN=[0 4];
JSPAN=[0 5];

%rule for jumps
% rule = 1 priority for jumps
% rule = 2 priority for flow
% rule = 3 random priority
rule = 1;

options=odeset('RelTol', 1e-4, 'MaxStep',0.01);

[t j x] = HyEQsolver(@f,@g,@C,@D, x0, TSPAN, JSPAN, rule, options);

figure
subplot(4,1,1)
plotflows(t,j,x(:,1))
ylabel('Z')
subplot(4,1,2)
plotflows(t,j,x(:,2))
ylabel('Tau')
subplot(4,1,3)
plotflows(t,j,x(:,3))
ylabel('q')
subplot(4,1,4)
plotflows(t,j,x(:,4))
ylabel('Phi')

PLL with phase-frequency detector

The idea is to model a classic PLL with a Phase-Frequency Detector (PFD) as a hybrid system. This post documents the first attempt…

The diagram of the PLL is:

pll_1

The loop filter and the VCO can be almost perfectly modeled as linear subsystems (ignoring non-linearities due to circuit design, saturations, etc.).

The Loop Filter proposed for this analysis has a current input and a voltage output:

pll_2

The VCO has infinite input impedance. The filter’s equations are (after some manipulations):

pll_eq_1

where v_{1} and v_{2} are the capacitor voltages, relative to VCC/2. The filter’s output is v_{0}=v_{1}.

The VCO’s equation is:

pll_eq_2

where k_{V} is the control input to frequency output gain, and f_{0} is the free-running frequency (which is the output when the control voltage is VCC/2).

The reference input, \phi_{REF}, is modeled as an oscillator with fixed frequency f_{REF}:

pll_eq_3

The PFD as the phase detector is inherently hybrid. Its internal circuit is:

pll_3

Block \tau_{D} is a delay element (used to eliminate a dead-band in the phase comparator to reduce output jitter).

The following hypotheses are made to simplify the analysis:

  • the two current sources are considered to be perfectly matched.
  • the delays of the two flip-flops (especially the delays in the reset circuits) are considered to be perfectly matched.
  • the switches have zero ON resistance.
  • the real delays in the OR gate and in the FFs’ reset circuits are included in \tau_{D}.

As an example of operation, the PFD’s output is shown for a particular case:

pll_4

If we work with the phases modulo-2\pi, we can draw its state machine as follows:

pll_5

so, in a compact form,

pll_eq_4

To model the delay, we use a timer \tau, such that \dot{\tau}=1 when u=d=1.

Summarizing, the state vector will be:

pll_eq_5

The flow map f(x) will be:

pll_eq_6

The jump map g(x) will be:

pll_eq_7

where

pll_eq_8

Looking at the PFD’s state machine we can obtain the jump set, and, as a consequence, also the flow set:

pll_eq_9

As an example, we can use this model to simulate a circuit based on the 74HC9046 PLL. In this example, we use the following parameters:

pll_eq_10

We can design a stabilizing loop filter (C_{1}C_{2}R) based on the linearized model of the PFD. Looking at the PFD’s signals, and assuming frequency tracking, we can see that the linearized model of the PFD takes the form:

pll_eq_11

The loop gain will be:

pll_eq_12

where Z\left(s\right) is the filter’s impedance to ground (seen from its input). If we take:

pll_eq_13

we can see that the linearized model will be stable:

pll_6

The step response of the linearized system is:

pll_7

To simulate this system, we can set an initial condition such as:

pll_eq_14

This means that the initial state is:

  • capacitors charged to VCC/2
  • PFD’s mode is u=d=0.
  • the timer is reset to 0.
  • the initial reference phase is \pi.
  • the initial VCO phase is 0.

This implies that the VCO’s frequency will be de-synchronized from the reference. We expect to see a traditional non-linear PLL acquisition behavior.

The following code simulates the system, using HyEQ Lite.

run.m:

%condiciones iniciales:
%x = v1, v2, u, d, tau, fi_ref, fi_vco
x0 = [0; 0;  0; 0; 0;   pi;     0];

%horizonte de simulación
TSPAN = [0 10];
JSPAN = [0 2500];

% rule for jumps
% rule = 1 -> priority for jumps
% rule = 2 -> priority for flows
% rule = 3 -> random priority
rule = 1;

options = odeset('RelTol',1e-6,'MaxStep',.1);

%simular
[t j x] = HyEQsolver(@f, @g, @C, @D, x0, TSPAN, JSPAN, rule, options);

%graficar la frecuencia
Kv = 1E6;
f0 = 0.9E6;
figure
plot(t * 1E3,(Kv*x(:,1) + f0) / 1E6);
xlabel 't (ms)'
ylabel 'f (MHz)'

Flow map f.m:

function xdot = f(x)

    %parámetros: C1, C2, R, ip, Kv, fref
    ip = 1E-3;      %1mA
    C1 = 390E-9;    %390nF
    C2 = 3.9e-6;    %3.9uF
    R = 50;         %50ohms
    Kv = 1E6;
    fref = 1E6;     %1MHz
    f0 = 0.9E6;

    %x = v1, v2, u, d, tau, fi_ref, fi_vco
    v1 = x(1);
    v2 = x(2);
    u = x(3);
    d = x(4);
    tau = x(5);
    fi_ref = x(6);
    fi_vco = x(7);

    xdot =  [1/C1*(u-d)*ip - (v1-v2)/(R*C1);   (v1 - v2)/(R*C2);   0; 0; u*d; 2*pi*fref; 2*pi*(Kv*v1+f0)];

end

Jump map g.m:

function xplus = g(x)

    %parámetros: Tdelay
    Tdelay = 10E-9;  %10ns

    %x = v1, v2, u, d, tau, fi_ref, fi_vco
    v1 = x(1);
    v2 = x(2);
    u = x(3);
    d = x(4);
    tau = x(5);
    fi_ref = x(6);
    fi_vco = x(7);

    %calculo u, d, tau
    if (tau >= Tdelay) %asumo u=d=1
        u = 0;
        d = 0;
        tau = 0;
    elseif (u == 0) && (d == 0)
        if (fi_ref >= 2*pi)
            u = 1;
            d = 0;
        else %asumo fi_vco >= 2*pi
            u = 0;
            d = 1;
        end
    else
        u = 1;
        d = 1;
    end

    %calculo fi_*
    if ((fi_ref >= 2*pi) || (fi_vco >= 2*pi))
        fi_ref = fi_ref - 2*pi;
        fi_vco = fi_vco - 2*pi;
    end

    xplus = [v1; v2; u; d; tau; fi_ref; fi_vco];

end

Flow set C.m:

function b = C(x)
    b = 1 - D(x);
end

Jump set D.m:

function b = D(x)
    %parámetros: Tdelay
    Tdelay = 10E-9;  %10ns

    %x = v1, v2, u, d, tau, fi_ref, fi_vco
    v1 = x(1);
    v2 = x(2);
    u = x(3);
    d = x(4);
    tau = x(5);
    fi_ref = x(6);
    fi_vco = x(7);

    if (u == 0) && (d == 0) && ((fi_ref >= 2*pi) || (fi_vco >= 2*pi))
        b = 1;
    elseif (u == 0) && (d == 1) && (fi_ref >= 0)
        b = 1;
    elseif (u == 1) && (d == 0) && (fi_vco >= 0)
        b = 1;
    elseif (u == 1) && (d == 1) && (tau >= Tdelay)
        b = 1;
    else
        b = 0;
    end

end

We can plot some results (for example, the VCO’s output frequency):

pll_8

We can see that roughly until 0.5ms the PLL exhibits a non-linear behavior, and after that, the linearized model seems to be valid.

A zoom of the last part of the curve shows the switching:

pll_9

A zoom of the acquisition phase shows the non-linear behavior:

pll_10