風速時程模擬自回歸法空間20個點-AR模型

1、%風速時程模擬自回歸法空間20個點-AR模型%自回歸模型階p4，模擬空間20個點，時間步長ti0.1，頻率步長f0.001，%空間相干系數(shù)采用與頻率無關的shiotani相關系數(shù),脈動風速譜為Davenport譜 cleartick=0.005;v10=25; n=0.001:0.001:10;xn=1200*n./v10;s1=4*k*252*xn.2./n./(1+xn.2).(4/3); %Davenport譜%產(chǎn)生空間點坐標for i=1:20 x(i)=5+i； z(i)=8+i;end %求R矩陣syms fR0=zeros(20);for i=1:20 for j=i:20 H0

2、=inline(4*12002*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); endend R1=zeros(20);for i=1:20 for j=i:20 H1=inline(4*12002*f*k).*exp(-sqrt(dx2/502+dz2/602).*cos(2*pi*f*ti)./(1+(1200*f/v10).2).(4/3),f,k,dx,dz,ti,v10); k=0.

3、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); endend R2=zeros(20);for i=1:20 for j=i:20 H2=inline(4*12002*f*k).*exp(-sqrt(dx2/502+dz2/602).*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=

4、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); endend R3=zeros(20);for i=1:20 for j=i:20 H3=inline(4*12002*f*k).*exp(-sqrt(dx2/502+dz2/602).*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(

5、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); endend R4=zeros(20);for i=1:20 for j=i:20 H4=inline(4*12002*f*k).*exp(-sqrt(dx2/502+dz2/602).*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(

6、H4,0.001,10,0.001,0,k,dx,dz,ti,v10); R4(j,i)=R4(i,j); endend 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=AB; %此式相當于A*X=B，X(8020矩陣)為自回歸系數(shù) q1=X(1:20,:); %取X的第一個2020矩陣q2=X(20+1:2*20,:); %取X的第二個2020矩陣q3=X(2*20+1

7、:3*20,:); %取X的第三個2020矩陣q4=X(3*20+1:4*20,:); %取X的第四個2020矩陣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); %產(chǎn)生均值0，方差1的正態(tài)隨機數(shù)矩陣 endend V(1:20,1)=L*a(:,1);V(1:20,2)=-Q1*V(1:20,1)+L*a(:,

8、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 V1=V(1,:); %取第一點的風速 t=(1:2048)*0.1;figuresubplot(2,1,1);plot(t,V1,b-);xlabel(t(s);ylabel(v(t);axis(0 120 -40 40); %與目標譜進行比較power,freq=psd(V(1,