风速时程模拟自回归法空间20个点-AR模型

%风速时程模拟自回归法空间20个点-AR模型

%自回归模型阶p＝4，模拟空间20个点，时间步长ti＝0.1，频率步长f＝0.001，

%空间相干系数采用与频率无关的shiotani相关系数,脉动风速谱为Davenport谱

clear

tic

k=0.005;

v10=25;

n=0.001:0.001:10;

xn=1200*n./v10;

s1=4*k*25^2*xn.^2./n./(1+xn.^2).^(4/3); %Davenport谱

%产生空间点坐标

for i=1:20

x(i)=5+i；

z(i)=8+i;

end

%求R矩阵

syms f

R0=zeros(20);

for i=1:20

for j=i:20

H0=inline('(4*1200^2*f*k)./(1+(1200*f/v10).^2).^(4/3)','f','k','v10');

k=0.005; %地面粗糙度长度

v10=25;

R0(i,j)=quadl(H0,0.001,10,0.001,0,k,v10);

R0(j,i)=R0(i,j);

end

end

R1=zeros(20);

for i=1:20

for j=i:20

H1=inline('(4*1200^2*f*k).*exp(-sqrt(dx^2/50^2+dz^2/60^2)).*cos(2*pi*f*

ti)./(1+(1200*f/v10).^2).^(4/3)','f','k','dx','dz','ti','v10');

k=0.005;

ti=0.1; %时间步长

v10=25;

dx=x(i)-x(j);

dz=z(i)-z(j);

R1(i,j)=quadl(H1,0.001,10,0.001,0,k,dx,dz,ti,v10);

R1(j,i)=R1(i,j);

end

end

R2=zeros(20);

for i=1:20

for j=i:20

H2=inline('(4*1200^2*f*k).*exp(-sqrt(dx^2/50^2+dz^2/60^2)).*cos(2*pi*f*

2*ti)./(1+(1200*f/v10).^2).^(4/3)','f','k','dx','dz','ti','v10');

k=0.005;

ti=0.1;

v10=25;

dx=x(i)-x(j);

dz=z(i)-z(j);

R2(i,j)=quadl(H2,0.001,10,0.001,0,k,dx,dz,ti,v10);

R2(j,i)=R2(i,j);

end

end

R3=zeros(20);

for i=1:20

for j=i:20

H3=inline('(4*1200^2*f*k).*exp(-sqrt(dx^2/50^2+dz^2/60^2)).*cos(2*pi*f*

3*ti)./(1+(1200*f/v10).^2).^(4/3)','f','k','dx','dz','ti','v10');

k=0.005;

ti=0.1;

v10=25;

dx=x(i)-x(j);

dz=z(i)-z(j);

R3(i,j)=quadl(H3,0.001,10,0.001,0,k,dx,dz,ti,v10);

R3(j,i)=R3(i,j);

end

end

R4=zeros(20);

for i=1:20

for j=i:20

H4=inline('(4*1200^2*f*k).*exp(-sqrt(dx^2/50^2+dz^2/60^2)).*cos(2*pi*f*

4*ti)./(1+(1200*f/v10).^2).^(4/3)','f','k','dx','dz','ti','v10');

k=0.005;

ti=0.1;

v10=25;

dx=x(i)-x(j);

dz=z(i)-z(j);

R4(i,j)=quadl(H4,0.001,10,0.001,0,k,dx,dz,ti,v10);

R4(j,i)=R4(i,j);

end

end

Q1=zeros(20);

Q2=zeros(20);

Q3=zeros(20);

Q4=zeros(20);

A=[R0 R1 R2 R3;R1 R2 R3 R0;R2 R3 R0 R1;R3 R0 R1 R2]; %80X80矩阵

B=[R1;R2;R3;R4]; %80X20矩阵

X=A\B; %此式相当于A*X=B，X(80×20矩阵)为自回归系数Ψ

q1=X(1:20,:); %取X的第一个20×20矩阵

q2=X(20+1:2*20,:); %取X的第二个20×20矩阵

q3=X(2*20+1:3*20,:); %取X的第三个20×20矩阵

q4=X(3*20+1:4*20,:); %取X的第四个20×20矩阵

Q1=q1';

Q2=q2';

Q3=q3';

Q4=q4';

RN=R0+Q1*R1+Q2*R2+Q3*R3+Q4*R4;

%对RN 进行cholesky分解

L=zeros(20);

L=chol(RN);

L=L';

a=zeros(20,2048);

for i=1:20

for j=1:2048

a(i,j)=normrnd(0,1,1,1); %产生均值0，方差1的正态随机数矩阵end

end

V(1:20,1)=L*a(:,1);

V(1:20,2)=-Q1*V(1:20,1)+L*a(:,2);

V(1:20,3)=-(Q1*V(1:20,2)+Q2*V(1:20,1))+L*a(:,3);

V(1:20,4)=-(Q1*V(1:20,3)+Q2*V(1:20,2)+Q3*V(1:20,1))+L*a(:,4);

for t=5:2048

V(1:20,t)=-(Q1*V(1:20,t-1)+Q2*V(1:20,t-2)+Q3*V(1:20,t-3)+Q4*V(1:20,t-4))+L*a(:,t); end