DCNLProjectHeartbeatModel/lyapunov.m at main · AndreaMazzera/DCNLProjectHeartbeatModel · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
%% Analisi di Stabilità
close all;

[xi1_mesh_Vdot1,xi2_mesh_Vdot1,xi1_meshV1,xi2_meshV1,V1,Vdot1,E1]=plot_lyapu(x_equ_diast,eps,T,xd,xs,0);
[xi1_mesh_Vdot2,xi2_mesh_Vdot2,xi1_meshV2,xi2_meshV2,V2,Vdot2,E2]=plot_lyapu(x_equ_sist,eps,T,xd,xs,1);

%% Plot derivate funzioni di Lyapunov
figure(2)
subplot(1,2,1)
hold on
grid on
box on
set(gca, 'FontSize', 14);
set(gcf, 'color', 'white', 'Position', [250 70 1100 500]);
view([-46 30])

plot3(x_equ_diast(1),x_equ_diast(2),0,'o','MarkerFaceColor','red','MarkerSize',15)
surf(xi1_mesh_Vdot1, xi2_mesh_Vdot1, Vdot1);
colorbar;
line([E1(1,1) E1(end,1)],[E1(1,2) E1(end,2)],'Color','k','LineWidth',5);
title('Stato Diastolico - $\dot{V}$','fontsize',14,'interpreter','latex')
xlabel('$x_1$', 'interpreter', 'latex');
ylabel('$x_2$', 'interpreter', 'latex');
zlabel('$\dot{V}$', 'interpreter', 'latex');

subplot(1,2,2)
hold on
grid on
box on
set(gca, 'FontSize', 14);
set(gcf, 'color', 'white');
view([-46 30])

plot3(x_equ_sist(1),x_equ_sist(2),0,'o','MarkerFaceColor','red','MarkerSize',15)
surf(xi1_mesh_Vdot2, xi2_mesh_Vdot2, Vdot2);
colorbar;
line([E2(1,1) E2(end,1)],[E2(1,2) E2(end,2)],'Color','k','LineWidth',5);
title('Stato Sistolico - $\dot{V}$','fontsize',14,'interpreter','latex')
xlabel('$x_1$', 'interpreter', 'latex');
ylabel('$x_2$', 'interpreter', 'latex');
zlabel('$\dot{V}$', 'interpreter', 'latex');
sgtitle('$\dot{V} = (x-x^*)^TP\dot{x}$', 'interpreter', 'latex')

%% Plot funzioni di Lyapunov
figure(3)
subplot(1,2,1)
hold on
grid on
box on
set(gca, 'FontSize', 14);
set(gcf, 'color', 'white', 'Position', [250 70 1100 500]);
view([-50 11])

plot3(x_equ_diast(1),x_equ_diast(2),0,'o','MarkerFaceColor','red','MarkerSize',15)
meshc(xi1_meshV1,xi2_meshV1,V1)
title('Stato Diastolico - $V$','fontsize',14,'interpreter','latex')
xlabel('$x_1$', 'interpreter', 'latex');
ylabel('$x_2$', 'interpreter', 'latex');
zlabel('$V$', 'interpreter', 'latex');

subplot(1,2,2)
hold on
grid on
box on
set(gca, 'FontSize', 14);
set(gcf, 'color', 'white');
view([-50 11])

plot3(x_equ_sist(1),x_equ_sist(2),0,'o','MarkerFaceColor','red','MarkerSize',15)
meshc(xi1_meshV2,xi2_meshV2,V2)
title('Stato Sistolico - $V$','fontsize',14,'interpreter','latex')
xlabel('$x_1$', 'interpreter', 'latex');
ylabel('$x_2$', 'interpreter', 'latex');
zlabel('$V$', 'interpreter', 'latex');

sgtitle('$V = \frac{1}{2}(x-x^*)^TP(x-x^*)$', 'interpreter', 'latex')


function [xi1_mesh_Vdot,xi2_mesh_Vdot,xi1_mesh_V,xi2_mesh_V,V,Vdot,E]=plot_lyapu(xequ,eps,T,xd,xs,u)

% Funzione che calcolo la meshgrid della V e Vdot, i quali di questi ultime
% infine calcola la regione E in cui la Vdot si annulla (retta nera)

    side = 0.5;
    step = 0.05;

    %P=eye(2,2);
    P=diag([0.1 0.03]);

    xi1_vec = xequ(1) - side : step : xequ(1) + side;
    xi2_vec = xequ(2) - side : step : xequ(2) + side;

    n1 = length(xi1_vec);
    n2 = length(xi2_vec);

    Vdot = zeros(n1, n2);
    V = zeros(n1, n2);

    for i = 1 : n1
        for j = 1 : n2
            xi=[xi1_vec(i); xi2_vec(j)];
            xi_dot=heartbeat_fun(xi,eps,T,xd,xs,u);
            Vdot(j,i) = Vdot_fun(xi, xequ, xi_dot, P);
            V(j,i)=0.5*(xi-xequ)'*P*(xi-xequ);
         end
    end

    a=find(Vdot >=0);
    size(a);

    [xi1_mesh_Vdot, xi2_mesh_Vdot] = meshgrid(xi1_vec, xi2_vec);
    s = surf(xi1_mesh_Vdot, xi2_mesh_Vdot, Vdot);

    E=zeros(length(s.XData),2);
    for i=1:length(s.XData)
        for j=1:length(s.YData)
            if s.ZData(i,j)==0
                E(j,1)=s.XData(i,j);
                E(j,2)=s.YData(i,j);
            end
        end
    end
    [xi1_mesh_V, xi2_mesh_V] = meshgrid(xi1_vec, xi2_vec);

% Retta x1=1.24
    close all
end