《量化投资:以MATLAB为工具》

MATLAB技术论坛

 找回密码
 注册账号
查看: 919|回复: 0
收起左侧

[提问] 求助大神 给看看这个 程序

[复制链接]
发表于 2015-7-11 21:12:13 | 显示全部楼层 |阅读模式
clear
tic
x=[0 0 0 0 0 0 0 0 0 0 0 0 0 0 ];
x=x';
y=[4.02725 3.3011 2.8372 2.4234 2.06155 1.8402 1.7439 1.645 1.546 1.44 1.3349 1.2369 1.1397 1.0046];
y=y';
z=[3.75 9.75 13.75 17.25 21.5 25.25 27.825 30.475 33.125 35.95 38.775 41.4 44 46.05];
z=z';
h=14;%节点个数
tqd=1;
v10=29.665;%10米处风速
k=1:h;
v(k)=v10*(z(k)/10).^0.1;  %任意高度处的风速
Tstep=6000;%步长的个数
dt=0.1;%步长的大小
f=0.01:0.01:10;%脉动风频率
u=0.005*v10^2;
P=4;
t=(1:6000)*dt;
nt=length(t);

%以下是求RO~R4
R0=zeros(h);  
for i=1:h
    for j=i:h
        dx=x(i)-x(j);
        dy=y(i)-y(j);
        dz=z(i)-z(j);
        h1=z(i)*f/v10/(z(i)/10)^0.16;
        h2=z(j)*f/v10/(z(j)/10)^0.16;
        w1=f.*(1+50*h1).^(5/3);
        w2=f.*(1+50*h2).^(5/3);                          
        s1=200*h1*u./w1;  %求kaimal谱
        s2=200*h2*u./w2;
        xiangguanhanshu=exp(-sqrt(dx.^2./50^2+dy.^2./50^2+dz.^2./60^2)).*cos(2*pi.*f*0.*dt);   %求相关函数
        y=sqrt(s1.*s2).*xiangguanhanshu;  %求互功率谱
        R0(i,j)=trapz(f,y);
        R0(j,i)=R0(i,j);
    end
end

R1=zeros(h);
for i=1:h
    for j=i:h
        dx=x(i)-x(j);
        dy=y(i)-y(j);
        dz=z(i)-z(j);
        h1=z(i)*f/v10/(z(i)/10)^0.16;
        h2=z(j)*f/v10/(z(j)/10)^0.16;
        w1=f.*(1+50*h1).^(5/3);
        w2=f.*(1+50*h2).^(5/3);
        s1=200*h1*u./w1;
        s2=200*h2*u./w2;
        xiangguanhanshu=exp(-sqrt(dx.^2./50^2+dy.^2./50^2+dz.^2./60^2)).*cos(2*pi.*f*1.*dt);
        y=sqrt(s1.*s2).*xiangguanhanshu;
        R1(i,j)=trapz(f,y);
        R1(j,i)=R1(i,j);
    end
end
   
R2=zeros(h);
for i=1:h
    for j=i:h
        dx=x(i)-x(j);
        dy=y(i)-y(j);
        dz=z(i)-z(j);
        h1=z(i)*f/v10/(z(i)/10)^0.16;
        h2=z(j)*f/v10/(z(j)/10)^0.16;
        w1=f.*(1+50*h1).^(5/3);
        w2=f.*(1+50*h2).^(5/3);
        s1=200*h1*u./w1;
        s2=200*h2*u./w2;
        xiangguanhanshu=exp(-sqrt(dx.^2./50^2+dy.^2./50^2+dz.^2./60^2)).*cos(2*pi.*f*2.*dt);
        y=sqrt(s1.*s2).*xiangguanhanshu;
        R2(i,j)=trapz(f,y);
        R2(j,i)=R2(i,j);
    end
end

R3=zeros(h);
for i=1:h
    for j=i:h
        dx=x(i)-x(j);
        dy=y(i)-y(j);
        dz=z(i)-z(j);
        h1=z(i)*f/v10/(z(i)/10)^0.16;
        h2=z(j)*f/v10/(z(j)/10)^0.16;
        w1=f.*(1+50*h1).^(5/3);
        w2=f.*(1+50*h2).^(5/3);
        s1=200*h1*u./w1;
        s2=200*h2*u./w2;
        xiangguanhanshu=exp(-sqrt(dx.^2./50^2+dy.^2./50^2+dz.^2./60^2)).*cos(2*pi.*f*3.*dt);
        y=sqrt(s1.*s2).*xiangguanhanshu;
        R3(i,j)=trapz(f,y);
        R3(j,i)=R3(i,j);
    end
end

R4=zeros(h);
for i=1:h
    for j=i:h
        dx=x(i)-x(j);
        dy=y(i)-y(j);
        dz=z(i)-z(j);
        h1=z(i)*f/v10/(z(i)/10)^0.16;
        h2=z(j)*f/v10/(z(j)/10)^0.16;
        w1=f.*(1+50*h1).^(5/3);
        w2=f.*(1+50*h2).^(5/3);
        s1=200*h1*u./w1;
        s2=200*h2*u./w2;
        xiangguanhanshu=exp(-sqrt(dx.^2./50^2+dy.^2./50^2+dz.^2./60^2)).*cos(2*pi.*f*4.*dt);
        y=sqrt(s1.*s2).*xiangguanhanshu;
        R4(i,j)=trapz(f,y);
        R4(j,i)=R4(i,j);
    end
end

Q1=zeros(h); %14*14阶零阵
Q2=zeros(h);
Q3=zeros(h);
Q4=zeros(h);
A=[-R0 -R1 -R2 -R3;-R1 -R0 -R1 -R2;-R2 -R1 -R0 -R1;-R3 -R2 -R1 -R0];
B=[R1;R2;R3;R4];
X=inv(A)*B;
q1=X(1:h,:);
q2=X(h+1:2*h,:);
q3=X(2*h+1:3*h,:);
q4=X(3*h+1:4*h,:);
Q1=q1';
Q2=q2';
Q3=q3';
Q4=q4';
RN=R0+Q1*R1+Q2*R2+Q3*R3+Q4*R4;  %求RN,协方差矩阵


L=chol(RN); %对RN柯西分解,求L
L=L';
a=zeros(h,nt);
a=normrnd(0,1,h,nt);  %构造 以0为均值,以1为方差,h*nt阶彼此相互独立的 的正态分布的随机函数
v(1:h,1)=L*a(:,1);
v(1:h,2)=-Q1*v(1:h,1)+L*a(:,2);
v(1:h,3)=-(Q1*v(1:h,2)+Q2*v(1:h,1))+L*a(:,3);
v(1:h,4)=-(Q1*v(1:h,3)+Q2*v(1:h,2)+Q3*v(1:h,1))+L*a(:,4);
for it=5:nt  %14个空间点相关脉动风速时程列向量的AR模型
    v(1:h,it)=-(Q1*v(1:h,it-1)+Q2*v(1:h,it-2)+Q3*v(1:h,it-3)+Q4*v(1:h,it-4))+L*a(:,it);
end
v;
%以下是平均风加脉动风,得总风速
fs1=v(1,:)+v10*(z(1)/10).^0.16;
fs2=v(2,:)+v10*(z(2)/10).^0.16;
fs3=v(3,:)+v10*(z(3)/10).^0.16;
fs4=v(4,:)+v10*(z(4)/10).^0.16;
fs5=v(5,:)+v10*(z(5)/10).^0.16;
fs6=v(6,:)+v10*(z(6)/10).^0.16;
fs7=v(7,:)+v10*(z(7)/10).^0.16;
fs8=v(8,:)+v10*(z(8)/10).^0.16;
fs9=v(9,:)+v10*(z(9)/10).^0.16;
fs10=v(10,:)+v10*(z(10)/10).^0.16;
fs11=v(11,:)+v10*(z(11)/10).^0.16;
fs12=v(12,:)+v10*(z(12)/10).^0.16;
fs13=v(13,:)+v10*(z(13)/10).^0.16;
fs14=v(14,:)+v10*(z(14)/10).^0.16;
fs=[fs1' fs2' fs3' fs4' fs5' fs6' fs7' fs8' fs9' fs10' fs11' fs12' fs13' fs14'];


v1=v(1,:);  %1,3,5点脉动风时程曲线对比
v3=v(3,:);
v5=v(5,:);
figure(1);
plot(t,v1,'b-');%画出v1时程
xlabel('t/s');
ylabel('v/(m/s)');
axis([0 600 -30 30]);  %x轴的范围0~600,y轴的范围-30~30

varr=var(v1) %求v1的方差
sumn=sum(v1)/nt%求v1的均值
[power,freq]=psd(v1,6144,10,boxcar(6072),0,'mean');%做功率谱密度函数,v1是信号,6144是快速傅里叶变换点数,10是采样频率,boxca矩形窗函数,6072是长度,
power=power*2*0.1;
figure(2);
i=tqd;  %以高度1处的风谱为目标普
u=0.005*v10^2;%式3-10  0.05为地面粗糙度系数  
h1=z(i).*f/v10./(z(i)/10)^0.16;  %同R0的注释
w1=f.*(1+50*h1).^(5/3);
s1=200*h1*u./w1;  %tqd点的kaimal谱
sn=s1';
loglog(freq,power,'r-',f,sn,'g-'); %x轴y轴均为对数函数, 绘制模拟普与目标普
xlabel('f/HZ');
ylabel('功率谱/(m^2/s)');
legend('模拟普','目标普');



figure(3);
[crra,cca]=xcorr(v(26,1:nt),6000);
[crrc,ccc]=xcorr(v(30,1:nt),6000);
[crr1,cc1]=xcorr(v(30,1:nt),v(26,1:nt),6000);
[crrd,ccd]=xcorr(v(30,1:nt),6000);
[crrf,ccf]=xcorr(v(78,1:nt),6000);
[crr2,cc2]=xcorr(v(78,1:nt),v(30,1:nt),6000);
d1=sqrt(crra(6001).*crrc(6001));
d2=sqrt(crrd(6001).*crrf(6001));
plot(cc1*0.1,crr1/d1,'r:',cc1*0.1,crr1/d2,'g--')
xlabel('t/s');
ylabel('互相关函数');

figure(4);
vtqd1=v10*(z(1)/10).^0.16;
v1m=v1+vtqd1;
plot(t,v1m,'b-');
xlabel('t/s');
ylabel('v/(m/s)');
legend('1点风速')
axis([0 600 0 50]);
%toc


%以下是各段的风荷载标准值
vv=fs.^2;
w0=vv./1600;
w01=w0(:,1)*2.59489*6.13849;
w02=w0(:,2)*2.555908*3.62265;
w03=w0(:,3)*2.519016*2.789576;
w04=w0(:,4)*2.485348*2.679088;
w05=w0(:,5)*2.383264*4.19356;
w06=w0(:,6)*2.30582*6.29741;
w07=w0(:,7)*2.397334*1.86239;
w08=w0(:,8)*2.381338*1.82649;
w09=w0(:,9)*2.37902*1.79619;
w010=w0(:,10)*2.32926*8.18112;
w011=w0(:,11)*2.425828*1.325464;
w012=w0(:,12)*2.401798*1.282566;
w013=w0(:,13)*2.375748*1.25946;
w014=w0(:,14)*2.28956*9.02759;
%w00=vpa(w00);%控制有效数字


%以下是各段的风荷载标准值平分到节点上
l01=w01./4;
l02=w02./4;
l03=w03./4;
l04=w04./4;
l05=w05./4;
l06=w06./4;
l07=w07./4;
l08=w08./4;
l09=w09./4;
l010=w010./4;
l011=w011./4;
l012=w012./4;
l013=w013./4;
l014=w014./4;

hz1=(l01+l02)./4;
hz2=(l01+l02)./4;
hz3=(l02+l03)./4;
hz4=(l02+l03)./4;
hz5=(l03+l04)./4;
hz6=(l03+l04)./4;
hz7=(l04+l05)./4;
hz8=(l04+l05)./4;
hz9=(l05+l06)./4;
hz10=(l05+l06)./4;
hz11=(l06+l07)./4;
hz12=(l06+l07)./4;
hz13=(l07+l08)./4;
hz14=(l07+l08)./4;
hz15=(l08+l09)./4;
hz16=(l08+l09)./4;
hz17=(l09+l010)./4;
hz18=(l09+l010)./4;
hz19=(l010+l011)./4;
hz20=(l010+l011)./4;
hz21=(l011+l012)./4;
hz22=(l011+l012)./4;
hz23=(l012+l013)./4;
hz24=(l012+l013)./4;
hz25=(l013+l014)./4;
hz26=(l013+l014)./4;
hz27=l014./4;
hz28=l014./4;
hz=1000*[hz1 hz2 hz3  hz4 hz5 hz6 hz7 hz8 hz9 hz10 hz11 hz12 hz13 hz14 hz15 hz16 hz17 hz18 hz19 hz20 hz21 hz22 hz23 hz24 hz25 hz26 hz27 hz28];
hz=hz(1:600,:);  %1分钟
ls=(1:28);
hh=1;
hz=[hz(1:hh-1,:);ls;hz(hh:end,:)];
hs=(0:600);
hz=[hs',hz];
save('hz30.txt', 'hz', '-ascii');%存为txt格式







??? Index exceeds matrix dimensions.运行结果出现这个了 不知道哪里的错误   求助大神帮忙
您需要登录后才可以回帖 登录 | 注册账号

本版积分规则

QQ|网站地图|MATLAB技术论坛|Simulink仿真论坛 ( 蜀ICP备19014457号-2 

GMT+8, 2022-8-18 19:02 , Processed in 0.052181 second(s), 22 queries , Gzip On.

Powered by Discuz! X3.4

© 2001-2017 Comsenz Inc.

快速回复 返回顶部 返回列表