# Sample-and-hold control

Consider the following sample-and-hold architecture:

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:

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

# 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:

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


# 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

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

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.

% 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

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

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

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.

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$.

# 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:

Entre 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:

Simulació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.

Figura 1: trayectoria de la variable $Z$ con respecto al tiempo híbrido.

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

Figura 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:

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:

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

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:

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}$:

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

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:

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

so, in a compact form,

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

Summarizing, the state vector will be:

The flow map f(x) will be:

The jump map g(x) will be:

where

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

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:

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:

The loop gain will be:

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

we can see that the linearized model will be stable:

The step response of the linearized system is:

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

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):

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:

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