clear all; close all; clc;
a1=0.1;a2=0.2;
omega1=5;omega2=4;
G=1;C12=0.01;C21=0.02;
dt=0.01;
x1(1)=0.5;
y1(1)=0.5;
x2(1)=1;
y2(1)=1;
for i=2:1000
x1(i)=x1(i-1)+((a1-x1(i-1)^2-y1(i-1)^2)*x1(i-1)-omega1*y1(i-1)+G*C12*(x2(i-1)-x1(i-1)))*dt;
y1(i)=y1(i-1)+((a1-x1(i-1)^2-y1(i-1)^2)*y1(i-1)+omega1*x1(i-1)+G*C12*(y2(i-1)-y1(i-1)))*dt;
x2(i)=x2(i-1)+((a2-x2(i-1)^2-y2(i-1)^2)*x2(i-1)-omega2*y2(i-1)+G*C21*(x1(i-1)-x2(i-1)))*dt;
y2(i)=y2(i-1)+((a2-x2(i-1)^2-y2(i-1)^2)*y2(i-1)+omega2*x2(i-1)+G*C21*(y1(i-1)-y2(i-1)))*dt;
end
N_measurements=31;
N_basis=31;
index=randi([100,999],1,N_measurements);
Xdot1=zeros([1,N_measurements]);
Ydot1=zeros([1,N_measurements]);
Xdot2=zeros([1,N_measurements]);
Ydot2=zeros([1,N_measurements]);
for ni=1:N_measurements
Xdot1(ni)=(x1(index(ni)+1)-x1(index(ni)))/dt;
Ydot1(ni)=(y1(index(ni)+1)-y1(index(ni)))/dt;
Xdot2(ni)=(x2(index(ni)+1)-x2(index(ni)))/dt;
Ydot2(ni)=(y2(index(ni)+1)-y2(index(ni)))/dt;
end
M=zeros([N_measurements,N_basis]);
for i=1:N_measurements
for j=1:N_basis
if j==1
M(i,j)=1;
elseif j==2
M(i,j)=x1(index(i));
elseif j==3
M(i,j)=y1(index(i));
elseif j==4
M(i,j)=x2(index(i));
elseif j==5
M(i,j)=y2(index(i));
elseif j==6
M(i,j)=x1(index(i))^2;
elseif j==7
M(i,j)=x2(index(i))^2;
elseif j==8
M(i,j)=y1(index(i))^2;
elseif j==9
M(i,j)=y2(index(i))^2;
elseif j==10
M(i,j)=x1(index(i))*x2(index(i));
elseif j==11
M(i,j)=x1(index(i))*y1(index(i));
elseif j==12
M(i,j)=x1(index(i))*y2(index(i));
elseif j==13
M(i,j)=x2(index(i))*y1(index(i));
elseif j==14
M(i,j)=x2(index(i))*y2(index(i));
elseif j==15
M(i,j)=y1(index(i))*y2(index(i));
elseif j==16
M(i,j)=x1(index(i))^3;
elseif j==17
M(i,j)=y1(index(i))^3;
elseif j==18
M(i,j)=x2(index(i))^3;
elseif j==19
M(i,j)=y2(index(i))^3;
elseif j==20
M(i,j)=x1(index(i))^2*x2(index(i));
elseif j==21
M(i,j)=x1(index(i))^2*y1(index(i));
elseif j==22
M(i,j)=x1(index(i))^2*y2(index(i));
elseif j==23
M(i,j)=x2(index(i))^2*x1(index(i));
elseif j==24
M(i,j)=x2(index(i))^2*y1(index(i));
elseif j==25
M(i,j)=x2(index(i))^2*y2(index(i));
elseif j==26
M(i,j)=y1(index(i))^2*x1(index(i));
elseif j==27
M(i,j)=y1(index(i))^2*x2(index(i));
elseif j==28
M(i,j)=y1(index(i))^2*y2(index(i));
elseif j==29
M(i,j)=y2(index(i))^2*x1(index(i));
elseif j==30
M(i,j)=y2(index(i))^2*x2(index(i));
else j==31
M(i,j)=y2(index(i))^2*y1(index(i));
end
end
end
Norms=zeros([1,N_basis]);
for j=1:N_basis
Norms(j)=norm(M(:,j));
end
for i=1:N_measurements
for j=1:N_basis
M(i,j)=M(i,j)/Norms(j);
end
end
lambda=0.001;
Xdot1=Xdot1';
cd 'E:\MATLAB\bin\cvx\cvx'
cvx_setup
cvx_begin
variable ksaix1(31);
minimize(norm(M*ksaix1-Xdot1)+ lambda*norm(ksaix1,1));
cvx_end
Ydot1=Ydot1';
cvx_setup
cvx_begin
variable ksaiy1(31);
minimize(norm(M*ksaiy1-Ydot1)+ lambda*norm(ksaiy1,1));
cvx_end
Xdot2=Xdot2';
cvx_setup
cvx_begin
variable ksaix2(31);
minimize(norm(M*ksaix2-Xdot2)+ lambda*norm(ksaix2,1));
cvx_end
Ydot2=Ydot2';
cvx_setup
cvx_begin
variable ksaiy2(31);
minimize(norm(M*ksaiy2-Ydot2)+ lambda*norm(ksaiy2,1));
cvx_end
for j=1:N_basis
ksaix1(j)=ksaix1(j)/Norms(j);
ksaiy1(j)=ksaiy1(j)/Norms(j);
ksaix2(j)=ksaix2(j)/Norms(j);
ksaiy2(j)=ksaiy2(j)/Norms(j);
end
ksaix1
ksaiy1
ksaix2
ksaiy2
figure
hold on
plot(x1,'r')
plot(x2,'g')
x1_rec(1)= x1(1);
x2_rec(1)= x2(1);
y1_rec(2)= y1(1);
y2_rec(2)= y2(1);
a1=ksaix1(1);
a2=ksaix1(2);
a3=ksaix1(3);
a4=ksaix1(4);
a5=ksaix1(5);
a6=ksaix1(6);
a7=ksaix1(7);
a8=ksaix1(8);
a9=ksaix1(9);
a10=ksaix1(10);
a11=ksaix1(11);
a12=ksaix1(12);
a13=ksaix1(13);
a14=ksaix1(14);
a15=ksaix1(15);
a16=ksaix1(16);
a17=ksaix1(17);
a18=ksaix1(18);
a19=ksaix1(19);
a20=ksaix1(20);
a21=ksaix1(21);
a22=ksaix1(22);
a23=ksaix1(23);
a24=ksaix1(24);
a25=ksaix1(25);
a26=ksaix1(26);
a27=ksaix1(27);
a28=ksaix1(28);
a29=ksaix1(29);
a30=ksaix1(30);
a31=ksaix1(31);
b1=ksaiy1(1);
b2=ksaiy1(2);
b3=ksaiy1(3);
b4=ksaiy1(4);
b5=ksaiy1(5);
b6=ksaiy1(6);
b7=ksaiy1(7);
b8=ksaiy1(8);
b9=ksaiy1(9);
b10=ksaiy1(10);
b11=ksaiy1(11);
b12=ksaiy1(12);
b13=ksaiy1(13);
b14=ksaiy1(14);
b15=ksaiy1(15);
b16=ksaiy1(16);
b17=ksaiy1(17);
b18=ksaiy1(18);
b19=ksaiy1(19);
b20=ksaiy1(20);
b21=ksaiy1(21);
b22=ksaiy1(22);
b23=ksaiy1(23);
b24=ksaiy1(24);
b25=ksaiy1(25);
b26=ksaiy1(26);
b27=ksaiy1(27);
b28=ksaiy1(28);
b29=ksaiy1(29);
b30=ksaiy1(30);
b31=ksaiy1(31);
c1=ksaix2(1);
c2=ksaix2(2);
c3=ksaix2(3);
c4=ksaix2(4);
c5=ksaix2(5);
c6=ksaix2(6);
c7=ksaix2(7);
c8=ksaix2(8);
c9=ksaix2(9);
c10=ksaix2(10);
c11=ksaix2(11);
c12=ksaix2(12);
c13=ksaix2(13);
c14=ksaix2(14);
c15=ksaix2(15);
c16=ksaix2(16);
c17=ksaix2(17);
c18=ksaix2(18);
c19=ksaix2(19);
c20=ksaix2(20);
c21=ksaix2(21);
c22=ksaix2(22);
c23=ksaix2(23);
c24=ksaix2(24);
c25=ksaix2(25);
c26=ksaix2(26);
c27=ksaix2(27);
c28=ksaix2(28);
c29=ksaix2(29);
c30=ksaix2(30);
c31=ksaix2(31);
d1=ksaiy2(1);
d2=ksaiy2(2);
d3=ksaiy2(3);
d4=ksaiy2(4);
d5=ksaiy2(5);
d6=ksaiy2(6);
d7=ksaiy2(7);
d8=ksaiy2(8);
d9=ksaiy2(9);
d10=ksaiy2(10);
d11=ksaiy2(11);
d12=ksaiy2(12);
d13=ksaiy2(13);
d14=ksaiy2(14);
d15=ksaiy2(15);
d16=ksaiy2(16);
d17=ksaiy2(17);
d18=ksaiy2(18);
d19=ksaiy2(19);
d20=ksaiy2(20);
d21=ksaiy2(21);
d22=ksaiy2(22);
d23=ksaiy2(23);
d24=ksaiy2(24);
d25=ksaiy2(25);
d26=ksaiy2(26);
d27=ksaiy2(27);
d28=ksaiy2(28);
d29=ksaiy2(29);
d30=ksaiy2(30);
d31=ksaiy2(31);
for i=2:32
y1_rec(i)= b1+b2*x1_rec(i-1)+b3*y1_rec(i-1)+b4*x2_rec(i-1)+b5*y2_rec(i-1)+b6*x1_rec(i-1)^2+b7*x2_rec(i-1)^2+b8*y1_rec(i-1)^2+b9*y2_rec(i-1)^2+b10*x1_rec(i-1)*x2_rec(i-1)+b11*x1_rec(i-1)*y1_rec(i-1)+b12*x1_rec(i-1)*y2_rec(i-1)+b13*x2_rec(i-1)*y1_rec(i-1)+b14*x2_rec(i-1)*y2_rec(i-1)+b15*y1_rec(i-1)*y2_rec(i-1)+b16*x1_rec(i-1)^3+b17*y1_rec(i-1)^3+b18*x2_rec(i-1)^3+b19*y2_rec(i-1)^3+b20*x1_rec(i-1)^2*x2_rec(i-1)+b21*x1_rec(i-1)^2*y1_rec(i-1)+b22*x1_rec(i-1)^2*y2_rec(i-1)+b23*x2_rec(i-1)^2*x1_rec(i-1)+b24*x2_rec(i-1)^2*y1_rec(i-1)+b25*x2_rec(i-1)^2*y2_rec(i-1)+b26*y1_rec(i-1)^2*x1_rec(i-1)+b27*y1_rec(i-1)^2*x2_rec(i-1)+b28*y1_rec(i-1)^2*y2_rec(i-1)+b29*y2_rec(i-1)^2*x1_rec(i-1)+b30*y2_rec(i-1)^2*x2_rec(i-1)+b31*y2_rec(i-1)^2*y1_rec(i-1);
y2_rec(i)= d1+d2*x1_rec(i-1)+d3*y1_rec(i-1)+d4*x2_rec(i-1)+d5*y2_rec(i-1)+d6*x1_rec(i-1)^2+d7*x2_rec(i-1)^2+d8*y1_rec(i-1)^2+d9*y2_rec(i-1)^2+d10*x1_rec(i-1)*x2_rec(i-1)+d11*x1_rec(i-1)*y1_rec(i-1)+d12*x1_rec(i-1)*y2_rec(i-1)+d13*x2_rec(i-1)*y1_rec(i-1)+d14*x2_rec(i-1)*y2_rec(i-1)+d15*y1_rec(i-1)*y2_rec(i-1)+d16*x1_rec(i-1)^3+d17*y1_rec(i-1)^3+d18*x2_rec(i-1)^3+d19*y2_rec(i-1)^3+d20*x1_rec(i-1)^2*x2_rec(i-1)+d21*x1_rec(i-1)^2*y1_rec(i-1)+d22*x1_rec(i-1)^2*y2_rec(i-1)+d23*x2_rec(i-1)^2*x1_rec(i-1)+d24*x2_rec(i-1)^2*y1_rec(i-1)+d25*x2_rec(i-1)^2*y2_rec(i-1)+d26*y1_rec(i-1)^2*x1_rec(i-1)+d27*y1_rec(i-1)^2*x2_rec(i-1)+d28*y1_rec(i-1)^2*y2_rec(i-1)+d29*y2_rec(i-1)^2*x1_rec(i-1)+d30*y2_rec(i-1)^2*x2_rec(i-1)+d31*y2_rec(i-1)^2*y1_rec(i-1);
x1_rec(i)= a1+a2*x1_rec(i-1)+a3*y1_rec(i-1)+a4*x2_rec(i-1)+a5*y2_rec(i-1)+a6*x1_rec(i-1)^2+a7*x2_rec(i-1)^2+a8*y1_rec(i-1)^2+a9*y2_rec(i-1)^2+a10*x1_rec(i-1)*x2_rec(i-1)+a11*x1_rec(i-1)*y1_rec(i-1)+a12*x1_rec(i-1)*y2_rec(i-1)+a13*x2_rec(i-1)*y1_rec(i-1)+a14*x2_rec(i-1)*y2_rec(i-1)+a15*y1_rec(i-1)*y2_rec(i-1)+a16*x1_rec(i-1)^3+a17*y1_rec(i-1)^3+a18*x2_rec(i-1)^3+a19*y2_rec(i-1)^3+a20*x1_rec(i-1)^2*x2_rec(i-1)+a21*x1_rec(i-1)^2*y1_rec(i-1)+a22*x1_rec(i-1)^2*y2_rec(i-1)+a23*x2_rec(i-1)^2*x1_rec(i-1)+a24*x2_rec(i-1)^2*y1_rec(i-1)+a25*x2_rec(i-1)^2*y2_rec(i-1)+a26*y1_rec(i-1)^2*x1_rec(i-1)+a27*y1_rec(i-1)^2*x2_rec(i-1)+a28*y1_rec(i-1)^2*y2_rec(i-1)+a29*y2_rec(i-1)^2*x1_rec(i-1)+a30*y2_rec(i-1)^2*x2_rec(i-1)+a31*y2_rec(i-1)^2*y1_rec(i-1);
x2_rec(i)= c1+c2*x1_rec(i-1)+c3*y1_rec(i-1)+c4*x2_rec(i-1)+c5*y2_rec(i-1)+c6*x1_rec(i-1)^2+c7*x2_rec(i-1)^2+c8*y1_rec(i-1)^2+c9*y2_rec(i-1)^2+c10*x1_rec(i-1)*x2_rec(i-1)+c11*x1_rec(i-1)*y1_rec(i-1)+c12*x1_rec(i-1)*y2_rec(i-1)+c13*x2_rec(i-1)*y1_rec(i-1)+c14*x2_rec(i-1)*y2_rec(i-1)+c15*y1_rec(i-1)*y2_rec(i-1)+c16*x1_rec(i-1)^3+c17*y1_rec(i-1)^3+c18*x2_rec(i-1)^3+c19*y2_rec(i-1)^3+c20*x1_rec(i-1)^2*x2_rec(i-1)+c21*x1_rec(i-1)^2*y1_rec(i-1)+c22*x1_rec(i-1)^2*y2_rec(i-1)+c23*x2_rec(i-1)^2*x1_rec(i-1)+c24*x2_rec(i-1)^2*y1_rec(i-1)+c25*x2_rec(i-1)^2*y2_rec(i-1)+c26*y1_rec(i-1)^2*x1_rec(i-1)+c27*y1_rec(i-1)^2*x2_rec(i-1)+c28*y1_rec(i-1)^2*y2_rec(i-1)+c29*y2_rec(i-1)^2*x1_rec(i-1)+c30*y2_rec(i-1)^2*x2_rec(i-1)+c31*y2_rec(i-1)^2*y1_rec(i-1);
end
figure
hold on
plot(x1_rec,'Y')
plot(x2_rec,'b')