程序里的figure2图形输出不了应该怎么改？？？

clear;

clc;

tic;

c1=[7.35         7.78         7.64         8.36         8.23         8.55         8.23         7.29         6.93         6.75         6.69         6.58         6.66         6.97         6.82         7.21         7.54         8.18         8.56         8.37         8.81         8.56         7.93         7.77

];

k1=[2.23         2.26         1.95         2.27         2.15         1.94         2.01         2.12         2.28         2.29         2.16         2.29         2.28         2.09         2.22         1.96         2.07         2.27         2.22         2.28         1.96         1.91         2.24         2.27

];

c2=[8.84         8.66         8.26         8.02         7.65         7.70         7.40         7.12         6.94         6.86         7.05         6.92         7.02         7.38         7.55         7.69         7.98         8.32         8.78         9.21         9.18         9.04         9.32         9.11

];

k2=[2.01         1.98         1.89         2.04         2.21         2.10         1.88         2.15         1.91         2.01         2.18         1.96         2.04         2.07         1.93         2.21         2.19         2.20         1.91         2.09         2.05         2.11         1.93         2.12

];

c3=[7.72         7.83         7.43         7.14         7.31         7.46         7.80         7.40         7.08         7.23         6.82         6.92         7.49         7.75         7.61         7.89         7.40         7.38         7.70         7.99         8.05         8.18         8.36         8.01

];

k3=[2.09         1.96         1.90         1.97         1.92         2.08         2.08         2.08         1.91         2.03         1.93         1.94         2.06         1.98         2.01         1.98         2.03         1.93         1.92         1.98         1.91         2.08         1.95         2.00

];

n=20000;

for kkk=1:24;

c=[c1(kkk);c2(kkk);c3(kkk)];%给定值

k=[k1(kkk);k2(kkk);k3(kkk)];%给定值

nwf=length(c);%数组长度(即行数或列数中的较大值)

rho=[1 0.76 0.55

0.76 1 0.27

0.55 0.27 1];%3个机组的相关性

y=copularnd('gaussian',rho,n);%形成满足guassian 正态copula 的随机数

y1=normcdf(y,0,1);%求y的标准正态分布累积概率

wv=zeros(length(c),n);

for i=1:length(c)

wv(i,:)=wblinv(y1(:,i),c(i),k(i));%求y1的逆weibull累积分布，即竖的相同的那个值

end

v1(kkk,:)=wv(1,:);

v2(kkk,:)=wv(2,:);

v3(kkk,:)=wv(3,:);%这3步将正态分布和weibull相互关联，并且进行转换

end;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

tt=3;%选择研究的时段

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

figure(1);

plot(v1(tt,:),v2(tt,:),'+');

vin=4;

vr=12;

vout=20;

pr=2;

Nw=15;

wr=pr*Nw;

for i=1:24;

for j=1:20000;

if v1(i,j)

w1(i,j)=0;

elseif v1(i,j)>vout;

w1(i,j)=0;

elseif (v1(i,j)>vr)&(v1(i,j)<=vout);%大于vr并且小于vout的含义

w1(i,j)=wr;

else

w1(i,j)=(v1(i,j)-vin)*wr/(vr-vin);%w1不满足上面的条件之外的选择

end;

if v2(i,j)

w2(i,j)=0;

elseif v2(i,j)>vout;

w2(i,j)=0;

elseif (v2(i,j)>vr)&(v2(i,j)<=vout);

w2(i,j)=wr;

else

w2(i,j)=(v2(i,j)-vin)*wr/(vr-vin);

end;

if v3(i,j)

w3(i,j)=0;

elseif v3(i,j)>vout;

w3(i,j)=0;

elseif (v3(i,j)>vr)&(v3(i,j)<=vout);

w3(i,j)=wr;

else

w3(i,j)=(v3(i,j)-vin)*wr/(vr-vin);%w1.w2，w3这三个程序重复了3遍看一个就可以

end;

end;

end;

wsum=w1+w2+w3;%总的风电的功率；

[NB_Num ,Gen_Num, PD_NB_Num, Time_Num, Branch_Num, Cost_Subsection_Num, Reference_NB, Gen_UC, PD_System_UC, PD_NB_UC, Branch_UC]...

=load_case_UC;

load=PD_System_UC(tt,2);%负荷均值,同时为预测值

f=load+load*0.015*randn(1,n);%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%r = μ + δ*randn(m,n)，μ为均值,δ为标准差(取负荷均值的1.5%)，m,n为生成矩阵维数，为负荷实际值

w=wsum(tt,:);%总风电功率的tt数段的实际数值；

w_forest=mean(w);;%总风电功率的tt数段的预测数值；

delt_PD=f-load;

delt_W=w-w_forest;

a=delt_PD-delt_W;%tt时段 负荷-风电偏差的总和，综合随机变量

[b xa]=hist(a,10000);%分区间统计，这里分5000个区间，可根据你的数据适当调整

%  num=numel(a);%num是数据的总个数

%  c=cumsum(b/num);%计算累计分布

%  plot(xa,c,'k-');%画图

%  title ('旋转备用的综合分布概率');

%

%  Mmean=0.9;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%0.99要进行改动

%  B=abs(c-Mmean);

%  [ca,indexx]=sort(B);

%  c(indexx(1))

%  Ac1=[xa(indexx(1))];%横坐标的值

F = ksdensity(a,xa,'function','cdf');%对数列进行统计

figure (2);

plot(xa,F);%输出图形

xlabel ('旋转备用');

ylabel ('概率分布值');

Ac1=interp1(F,xa,0.9);%0.9就是置信度，当F为0.9时横坐标的值为AC1，也就是所需要的上旋转备用，可以根据需要修改

[NB_Num ,Gen_Num, PD_NB_Num, Time_Num, Branch_Num, Cost_Subsection_Num, Reference_NB, Gen_UC, PD_System_UC, PD_NB_UC, Branch_UC]...

=load_case_UC;

ngen = Gen_Num;

% eq groups are separated by states rather than time periods

Aeq = [1 1 1];

beq = load-w_forest;

% Researve Cap. Constraints

Ares = [1 1 1];

brest = - Ac1  + sum(Gen_UC(:,5)) ;%3台火电机组的最大出力

% set of constraints

A = [ Ares];%矩阵，相当于例题中的A

B = [ brest];%最大值，旋转备用约束

vlb = Gen_UC(:,6);%第6列为最小出力

vub = Gen_UC(:,5);%第5列为最大出力

gencost_c = 2.*diag(Gen_UC(:,2)'); %将发电成本二次系数写成方阵，Gen_UC(:,2)是第2列的书

gencost_b = Gen_UC(:,3)'; %列矩阵，相当于例题中的f

gencost_a = sum(Gen_UC(:,4));%第4列的数值之和

[Pi,fval] = quadprog(gencost_c,gencost_b,A,B,Aeq,beq,vlb,vub,[],[]);%见例题，方括号里为不存在的数

disp('有功功率约束满足情况')

Aeq*Pi-beq

disp('旋转备用约束满足情况')

Ares*Pi - brest%Pi-(P(i,max)-Ri)≤0

disp('发电费用')

Pi'*gencost_c*Pi + gencost_b*Pi + gencost_a

toc;