clear all;clc;
global p1 C e 
load('data.mat');
p=data(2,:);a=data(1,:);
p=p/5;p1=350,C=350,e=0.01
m=200,U=3.01;L1=110.8;
%船舶電力負荷數據仿真。。。。
v=size(a,2);
fprintf(' 船舶電力負荷數據仿真 ...\n');
    st = cputime;
    x=[];
    for n=1:v
        x(:,n)=rand(m,1);
    end
i=1:v;
t=max(x(:,i));
d=ones(m,v);
for n=1:v
    data=t(n);
d(:,n)=data.*d(:,n);
end
for sn=1:v
    fprintf('%d ... \n',sn);
    e=0.001;ee=0.0001;
    r=find((d(:,sn)-x(:,sn))<e);
    rr=size(r,1);
while rr<p(sn)-1
    e=e+ee;
    r=find((d(:,sn)-x(:,sn))<e);
    rr=size(r,1);
end
s=size(r,1)-p(sn)+1;
if s>1
    e=e/10;
   ee=ee/10;
r=find((d(:,sn)-x(:,sn))<e);
rr=size(x,1);
while rr<p(sn)-1
    e=e+ee;
    r=find((d(:,sn)-x(:,sn))<e);
    rr=size(r,1);
end
end
x(r,sn)=1;
end
re=find(x<1);
x(re)=0;X=[];
for n=1:v
    X(:,n)=x(:,n).*a(n);
end
%輸出隨機電力負荷值Y，平均電力負荷值EY、
y=zeros(size(X,1),1);
for n=1:v
    y=y+X(:,n);
end
y=y+L1;
sumEY=0;sumDY=0;
for n=1:m
   sumEY=y(n)+sumEY;
end
EY=sumEY/m;
for n=1:m
    sumDY=(y(n)-EY).^2+sumDY;
end
DY=sumDY/m;
L=U*sqrt(DY)+EY+L1;
fprintf('Execution time : %4.1f seconds\n',cputime - st);
fprintf('電力負荷平均值EY： %4.1f kw \n',EY)
fprintf('電力負荷方差DY： %4.1f kw \n',DY)
fprintf('電力負荷最大值maxY: %4.1f kw \n',max(y))
fprintf('電力負荷最小值minY: %4.1f kw \n',min(y))
fprintf('~~~~~~~~~~~~~~~~~~~~~~~~~~ \n')
fprintf('置信度a=0.999的最大負荷預測值： \n')
fprintf('L=U(1-a)*sqrt(DY)+EY+L1= %4.1f kw \n',L)
h1=figure;
k=1:m;
plot (k,y(k),'r')
hold on
EY=EY*ones(1,m);
plot(k,EY)
fprintf('~~~~~~~~~~~~~~~~~~~~~~~~ \n')
[nsv,beta,error1,error2,trnY,X]=powerRecursive2(y) ;