Export result from fmincon

2 次查看（过去 30 天）
Mi Tung
Mi Tung 2020-10-6
Hi everyone. I;m using fmincon for solving optimization and i got a problem with him. Below show the code I'm using.
In which I want to export results "C_i" but in the result I got C_i = []. And also I want to display figure(1) and figure(2) in the results but I couldn't get it. Please help me for these. Thank you very much
Best regards
function [xsol,fval,history,searchdir] = runfmincon222
% Set up shared variables with OUTFUN
history.x = [];
history.fval = [];
searchdir = [];
% call optimization
Lp=27.31; % ????? ?????-?????????
Bp=7.08; % ?????? ?????-?????????
Tp=2.6; % ?????? ?????-?????????
Hp=3.42; % ?????? ?????-?????????
A=[];
b=[];
Aeq=[Bp -Lp 0 0 0 0 0 0 0; 0 Tp -Bp 0 0 0 0 0 0; 0 0 Hp -Tp 0 0 0 0 0; 0 0 0 0 0 -1.03 1 0 0];
beq=[0; 0; 0; 0];
lb=[24 5 2 3 0.52 0.76 0.82 60 10];
ub=[50 12 4 5 0.65 0.86 0.92 90 12];
x0=[27.12 6.86 2.67 3.47 0.56 0.83 0.86 77 11];
options = optimoptions(@fmincon,'OutputFcn',@outfun, 'TolCon',0, ...
'Display','iter','Algorithm','active-set', 'MaxFunctionEvaluations',100);
[xsol,fval] = fmincon('optim',x0,[],[],Aeq,beq,lb,ub,'nonlcond',options);
function stop = outfun(x,optimValues,state)
stop = false;
switch state
case 'init'
hold on
case 'iter'
% Concatenate current point and objective function
% value with history. x must be a row vector.
history.fval = [history.fval; optimValues.fval];
history.x = [history.x; x];
% Concatenate current search direction with
% searchdir.
searchdir = [searchdir;...
optimValues.searchdirection'];
% plot(x(1),x(2),'o');
% % Label points with iteration number and add title.
% % Add .15 to x(1) to separate label from plotted 'o'
% text(x(1)+.15,x(2),...
% num2str(optimValues.iteration));
% title('Sequence of Points Computed by fmincon');
case 'done'
hold off
otherwise
end
end
function f=optim(x)
global Pgr1 Wgr h0 ls k Tk D1 sr_11
gko=0.561; % óäåëüíûé âåñ êîðïóñà ñ îáîðóäîâàíèåì ê ðàñ÷åòíîìó âîäîèçìåùåíèþ
qob=0.023; % Èçìåðèòåëü ìàññ ïðîìûñëîâîãî è òåõíîëîãè÷åñêîãî îáîðóäîâàíèÿ
qm=0.02; % Èçìåðèòåëü ìàññ ãëàâíûõ è âñïîìîãàòåëüíûõ ìåõàíèçì
qz=0.03; % Èçìåðèòåëü ìàññ çàïàñà âîäîèçìåùåíèÿ
qsn=0.415; % Èçìåðèòåëü ìàññ ñíàáæåíèÿ
qb=0.015; % Èçìåðèòåëü ìàññû âîäÿíîé áàëëàñòû ïðè âîçâðàùåíèè â ïîðòó
qg1=24*165*10^-6; % Óäåëüíûé ðàñõîä òîïëèâà íà 1 ë.ñ..ñóòêè íà ïåðåõîäàõ (ò/ë.ñ..ñóòêè)
qg2=0.6*qg1; % Óäåëüíûé ðàñõîä òîïëèâà íà 1 ë.ñ..ñóòêè íà ïðîìûñëå (ò/ë.ñ..ñóòêè)
am=0.03; %Êîýôôèöèåíò ó÷èòûâàþùèé çàïàñû ñìàçî÷íîãî ìàñëà äëÿ ÃÄ
i1=0.7; %Êîýôôèöèåíò èñïîëüçîâàíèÿ ìîùíîñòè ÃÄ íà ïðîìûñëå
i2=0.85;%Êîýôôèöèåíò ó÷èòûâàþùèé ïðîäîëüæèòåëüíîñòü ðàáîòû ÃÄ íà ïðîìûñëå
k1=1.1; % Êîýôôèöèåíò ó÷èòûâàþùèé ðàñõîäà òîïëèâà è ìàñëà äëÿ âñïîìîãàòåëüíûõ ìåõàíèçì
kv=0.95; %Êîýôôèöèåíò ó÷èòûâàþùèé ïîòåðè ñêîðîñòè ñóäíà íà ïåðåõîäàõ
% kpa=1.1;%Êîýôôèöèåíò ó÷èòûâàþùèé ïîòåðè âðåìåíè íà ïðîìûñëå
kpa=0.75;
% kpa=0.69;
psr=5.25;%Ñðåäíèéñóòî÷íûé óëîâ
% vs=11; %Ñêîðîñòü ñâîáîäíîãî õîäà ñóäíà, óçë.
% Pgr1=77;%Ãðóçîïîäúåìíîñòü ñóäíà, ò
RL=375; %Äàëüíîñòü ðàéîíà ïðîìûñëà, ìèëü.
tst=2;%Âðåìÿ ñòîÿíêè â ïîðòó, ñóòêè
tper=RL/(24*x(9)*kv);%Âðåìÿ ïåðåõîäà íà ðàéîí ïðîìûñëà, ñóòêè
tpro=kpa*x(8)/psr; %Âðåìÿ íà ïðîìûñëå, ñóòêè
% tpro=11;
tmz=2; %âðåìÿ ìîðñêîãî çàïàñà, ñóòêè
tob=1;%Âðåìÿ äëÿ òåõíè÷åñêîãî îáñëóæèâàíèÿ çà ðåéñ, ñóòêè.
C=404; % êîýôôèöèåíò Ñ ïî À.È. Ðàêîâó
Avt=tst+2*tper+tpro+tmz+tob; % Àâòîíîìíîñòü ïëàâàíèÿ ñóäíà.
nek=12; %Êîëè÷åñòâî ýêèïàæà
uek=0.1;%Ò/÷åë.
uprv=80;%Ë/÷åë..ñóòêè
upro=3;%Êã/÷åë..ñóòêè
me=nek*uek;%ìàññà ýêèïàæà, ò
mprv=uprv*(Avt-tst-tob)*nek/1000;%ìàññà ïðåñíîé âîäû, ò
mpro=upro*(Avt-tst-tob)*nek/1000;%ìàññà ïðîâèçèÿ, ò
msn=me+mprv+mpro;%Ìàññà ñíàáæåíèÿ ïðè âûõîäå íà ïðîìûñëå
msn1=0.4*msn;%Ìàññà ñíàáæåíèÿ ïðè íà÷àëüíîì âîçâðàùåíèè â ïîðòó
%%% Êîýôôèöèåíòû â óðàâíåíèè ìàññ ñóäíà
a1=1-(gko+qz+qob+qb);
% a2=(-1)*0.4*qsn
a2=0;
a3=(-1)*(qm*x(9).^(4.5)/C+k1*(1+am)*(qg1*(tper+tmz))*x(9).^(4.5)/C);
% a4=(-1)*Pgr1;
a4=(-1)*(x(8)+msn1);
Y1=[a1 a2 a3 a4];
Y2=roots(Y1);
Y3=real(Y2).^3;
D1=Y3(1); % Íàèáîëüøåå âîäîèçìåùåíèå ñóäíà
Pt=D1^(1/3)*(k1*(1+am)*(qg1*(2*tper+tmz)+qg2*i1*i2*tpro)*x(9).^(4.5)/C);
Pt1=D1^(1/3)*k1*(1+am)*(qg1*(tper+tmz))*x(9).^(4.5)/C;
Pm=D1^(1/3)*qm*x(9).^(4.5)/C;
%%%Õîäêîñòü
global skor1 N Rtr
om=x(1)*(x(3)+x(2)/2)*(0.55+1.52*x(5));
v0=11*0.5144;
for i=1:0.1:1.5
v=v0*i;
fr=v/(9.81*x(1))^0.5;
if fr>0.36;
fr=0.35;
else
fr=v/(9.81*x(1))^0.5;
end
Re=v*x(1)/16.1;
CF0=10^3*0.455*(log10(Re*10^7)).^(-2.58);
Ca=0.5;
Cap=0.25;
Pdv=(x(1)*x(2)*x(3)*x(5)*1.025)^(1/3)*11^4.5/554.7;
% çàâèñèìîñòü Cr îò L/B
if x(1)/x(2)>3.5 & x(1)/x(2)<4.75;
a1=x(1)/x(2);
else
a1=3.857;
end
y1=[0.21 0.23 0.25 0.27 0.29 0.31 0.34 0.36 0.37];
e1=[3.5 3.75 4 4.25 4.5 4.75];
f1=[1.2 1.17 1.15 1.11 1.05 0
1.48 1.46 1.35 1.31 1.27 1.2
1.89 1.86 1.78 1.75 1.68 1.58
2.4 2.35 2.32 2.24 2.15 2.03
2.92 2.89 2.81 2.73 2.62 2.5
3.55 3.5 3.42 3.31 3.2 3.02
5 4.8 4.66 4.51 4.19 3.82
6.78 6.45 6.08 5.72 5.38 4.92
7.81 7.51 7 6.57 6.12 5.76];
c1=interp2(e1,y1,f1,a1,fr);
% çàâèñèìîñòü Cr îò B/T
if x(2)/x(3)>2.3 & x(2)/x(3)<3.2;
a2=x(2)/x(3);
else
a2=2.722;
end
y2=[0.21 0.23 0.25 0.27 0.29 0.31 0.35 0.36 0.37];
e2=[2.3 2.5 2.7 2.9 3.1 3.2];
f2=[1.28 1.24 1.19 1.15 1.11 1.08
1.57 1.55 1.52 1.49 1.45 1.41
1.98 1.97 1.94 1.89 1.79 1.71
2.47 2.45 2.41 2.34 2.25 2.2
2.98 2.97 2.94 2.89 2.79 2.71
3.58 3.55 3.51 3.46 3.32 3.29
4.28 4.23 4.16 4.08 3.99 3.92
5.48 5.41 5.34 5.21 5.04 5
7.18 7.06 6.95 6.86 6.74 6.69];
c2=interp2(e2,y2,f2,a2,fr);
% Çàâèñèìîñòü Cr îò beta
a3=x(7);
y3=[0.21 0.23 0.25 0.27 0.29 0.31 0.34 0.36 0.37];
e3=[0.7 0.75 0.8 0.85 0.9 0.95];
f3=[1.1 1.12 1.14 1.17 1.26 1.34
1.46 1.49 1.53 1.55 1.59 1.6
1.85 1.89 1.92 1.96 2 2.08
2.16 2.27 2.35 2.42 2.58 2.78
2.8 2.84 2.87 2.97 3.38 3.58
3.19 3.29 3.48 3.52 4 4.5
4.08 4.48 4.54 4.8 5.28 5.5
5.4 5.78 5.92 6.1 6.76 7
6.2 6.65 6.83 7.27 7.88 8.24];
c3=interp2(e3,y3,f3,a3,fr);
% Çàâèñèìîñòü Cr îò xc
a4=0;
y4=[0.21 0.23 0.25 0.27 0.29 0.31 0.34 0.36 0.37];
e4=[-0.025 -0.02 -0.015 -0.01 -0.005 0 0.005 0.01];
f4=[1.15 1.15 1.15 1.15 1.15 1.15 1.15 1.15
1.5 1.52 1.54 1.56 1.58 1.61 1.65 1.7
1.76 1.79 1.84 1.89 1.96 2 2.05 2.1
2.05 2.11 2.2 2.36 2.48 2.58 2.68 2.78
2.57 2.67 2.77 2.87 2.99 3.08 3.19 3.38
2.95 3.08 3.31 3.46 3.62 3.8 4 4.08
3.86 4.06 4.4 4.62 4.95 5.1 5.42 5.63
5.04 5.41 5.79 6.03 6.34 6.58 6.85 7
6.1 6.47 6.75 7 7.32 7.5 7.71 7.82];
c4=interp2(e4,y4,f4,a4,fr);
% Çàâèñèìîñòü Cr îò êîýôôèöèåíò ô è ô0
if x(5)/x(7)>0.55 & x(5)/x(7)<0.7;
a5=x(5)/x(7);
else
a5=0.637;
end
y5=[0.21 0.23 0.25 0.27 0.29 0.31 0.34 0.36 0.37];
e5=[0.55 0.575 0.6 0.625 0.65 0.675 0.7];
f5=[1.01 1.04 1.1 1.25 1.47 1.74 2
1.37 1.47 1.56 1.73 1.9 2.06 2.3
1.73 1.83 1.94 2.04 2.37 2.55 2.81
1.92 2.05 2.3 2.65 2.98 3.3 3.78
2.21 2.55 2.9 3.45 3.98 4.77 5.58
2.65 3 3.5 4.12 5 6 7.08
4 4.18 4.7 5.42 6.38 7.48 8.9
5.71 5.81 6.1 6.55 7.43 8.5 9.75
6.79 6.82 7 7.56 8.35 9.28 10.35];
c5=interp2(e5,y5,f5,a5,fr);
a6=0.6;
c6=interp2(e5,y5,f5,a6,fr);
Cr=(c1.*c2.*c3.*c4.*c5)./(c6.^4);
C=CF0+Ca+Cap+Cr';
R=C.*1.025.*(v.^2)*om*0.5;
% îïðåäåëåíèå òðåáóåìîé ìîùíîñòè ïðè òðàâëåíèè
Rtr=R/10^3;
t=0.078;
w=0.078;
n1=0.55;
n2=0.96;
n3=0.99;
n=n1*n2*n3*(1-t)/(1-w);
Pdtr=v*Rtr/n;
skor1=v/0.5144;
load Marineengine.mat;
% load Marineengine1.mat
En = [engine.Mass;engine.Power; engine.Length ];
ep = En(2,:)./Pdtr;
% rp = find(ep > 1.05 & ep <1.2)
rp = find(ep > 1.05);
[n4,m4] = min(En(1,rp));
ken=engine(rp(m4(1,end)));
Mdv = En(1,rp(m4(1,end)));
N = En(2,rp(m4(1,end)));
l = En(3,rp(m4(1,end)));
kp2=N/Pdtr;
model=ken.Name;
% if kp2 > 1.05 && kp2 <1.2
if kp2 > 1.05
N1=N;
break
end
end
%%% Ñòàòüÿ íàãðóçêè ìàññ ñóäíà
global Pb
Pko=gko*D1;
Pob=qob*D1;
Pz=qz*D1;
Dpor=Pko+Pob+Pm+Pz;
Psn=msn;
Pb=qb*D1;
DW=Pt+Psn+Pb;
%%%Êîìïîíîâêà ñóäíà
global lr1 le lmo lph lax ltr
ne=12; %Êîëè÷åñòâî ýêèïàæà, ÷åë.
ne1=5;%Êîëè÷åñòâî ÷ëåíîâ ýêèïàæà, ïðîæèâàþùèõ â ðóáêå ïåðâîãî ÿðóñà
ne2=ne-ne1;%êîëè÷åñòâî ÷ëåíîâ ýêèïàæà, ïðîæèâàþùèõ â æèëîì îòñåêå ¹6
ke=6; %Ïëîùàäü íà 1 ÷ëåíà ýêèïàæà, ì2/÷åë.
kb1=1.25;%Êîýôôèöèåíò ó÷èòûâàþùèé äëèíó âñïîìîãàòåëüíîãî ïîìåùåíèÿ â ðóáêå
krsw=1.3;%Êîýôôèöèåíò ó÷èòûâàþùèé äëèíó ïîìåùåíèÿ óñòàíîâêè ñèñòåìû îõëàæäåíèÿ;
kb2=1.25;%Êîýôôèöèåíò ó÷èòûâàþùèé äëèíó âñïîìîãàòåëüíîãî ïîìåùåíèÿ æèëîì îòñåêå ¹6
lr1=ke*ne1*kb2/(0.7*x(2)); %Äëèíà ðóáêè ïåðâîãî ÿðóñà
le=ke*ne2*krsw*kb2/x(2); %Äëèíà æèëîãî ïîìåùåíèÿ
kmo=1.98; %Êîýôôèöèåíò äëÿ äëèíû Ì.Î.
lmo=kmo*l;%Äëèíà Ì.Î., ì
lm=1.07*x(1); %Íàèáîëüøàÿ äëèíà ñóäíà, ì
lph=4.58; %Äëèíà ôîðïèêà, ì
lax=5.95; %Äëèíà àõòåðïèêà, ì
ltr=lm-lph-lax-lmo-le;%Äëèíà òðþìà, ì
hw=1;%âûñîòà âòðîãî äíà, ì
Wgr=ltr*x(2)*(x(4)-hw)*0.88;
%%% îñòîé÷èâîñòü ñóäíà
%Ìåòàöåíòðè÷åñêàÿ âûñîòà
kz=0.75;
h0=2*x(2)*sqrt((1.017-0.023)*x(6)/(x(6)+x(5))*(1.06-0.05)*x(6)^2/12/x(5))-kz*x(4);
% Îðäèíàòû äèàãðàììû ñòàòè÷åñêîé îñòîé÷èâîñòè
teta=10;
te=[10 20 30 40 50 60 70 80 90];
Hte1=[0.050 0.387 0.840 1.279 1.365 1.056 0.583 0.210 0];
Hte2=[-0.036 -0.241 -0.556 -0.722 -0.513 0.026 0.603 0.935 1.000];
Hte3=[0.151 0.184 0.081 -0.069 -0.155 -0.135 -0.062 -0.010 0];
Hte4=[0.010 0.062 0.135 0.155 0.069 -0.081 -0.184 -0.151 0];
for i=1:90;
fte1(i)=interp1(te,Hte1,i);
fte2(i)=interp1(te,Hte2,i);
fte3(i)=interp1(te,Hte3,i);
fte4(i)=interp1(te,Hte4,i);
lst(i)=0.5*x(2)*(1-0.96*x(3)/x(4))*fte1(i)+0.64*(1-1.032*x(3)/x(4))*x(4)*fte2(i)+1/11.4*...
(x(6)*x(2))^2/x(5)/x(3)*fte3(i)+1/11.4*(x(6)*x(2))^2/x(5)/x(3)*((0.64*(1-1.032*x(3)/x(4))*x(4))/...
(0.5*x(2)*(1-0.96*x(3)/x(4))))^3*fte4(i)-(kz*x(4)-x(6)/(x(6)+x(5))*x(3))*sin(i*pi/180);
end
ls=lst(max(lst)==lst);
tkr=find((lst)==ls);
% êðèòåðèÿ ïîãîäû
for i=10:90;
fte1(i)=interp1(te,Hte1,i);
fte2(i)=interp1(te,Hte2,i);
fte3(i)=interp1(te,Hte3,i);
fte4(i)=interp1(te,Hte4,i);
lstd(i)=0.5*x(2)*(1-0.96*x(3)/x(4))*fte1(i)+0.64*(1-1.032*x(3)/x(4))*x(4)*fte2(i)+1/11.4*...
(x(6)*x(2))^2/x(5)/x(3)*fte3(i)+1/11.4*(x(6)*x(2))^2/x(5)/x(3)*((0.64*(1-1.032*x(3)/x(4))*x(4))/...
(0.5*x(2)*(1-0.96*x(3)/x(4))))^3*fte4(i)-(kz*x(4)-x(6)/(x(6)+x(5))*x(3))*sin(i*pi/180);
end
ls0=0;
ls10=lstd(10);
ls20=lstd(20);
ls30=lstd(30);
ls40=lstd(40);
ls50=lstd(50);
ls60=lstd(60);
ls70=lstd(70);
ls80=lstd(80);
ls90=lstd(90);
X=0.0873;
l0=ls0*X;
l10=ls10*X;
l20=(2*ls10+ls20)*X;
l30=(2*ls10+2*ls20+ls30)*X;
l40=(2*ls10+2*ls20+2*ls30+ls40)*X;
l50=(2*ls10+2*ls20+2*ls30+2*ls40+ls50)*X;
l60=(2*ls10+2*ls20+2*ls30+2*ls40+2*ls50+ls60)*X;
l70=(2*ls10+2*ls20+2*ls30+2*ls40+2*ls50+2*ls60+ls70)*X;
l80=(2*ls10+2*ls20+2*ls30+2*ls40+2*ls50+2*ls60+2*ls70+ls80)*X;
l90=(2*ls10+2*ls20+2*ls30+2*ls40+2*ls50+2*ls60+2*ls70+2*ls80+ls90)*X;
kr=[0 10 20 30 40 50 60 70 80 90];
M1=[ls0 ls10 ls20 ls30 ls40 ls50 ls60 ls70 ls80 ls90];
M2=[l0 l10 l20 l30 l40 l50 l60 l70 l80 l90];
Fp0 = griddedInterpolant(kr,M1,'spline');
funp0 = @(t) Fp0(t);
Fp3=griddedInterpolant(kr,M2,'spline');
funp3 = @(t) Fp3(t);
Av=(0.96*1.15*x(1)*(1.68*x(4)-x(3))+2.4*8.2*0.65);
zp=2.99;
Pvtr=504;
lw1=Pvtr*zp*Av/(1000*9.81*x(1)*x(2)*x(3)*x(5)*1.025);
lw2=lw1*1.5;
ft1=[0 10 20 30];
ft2=[ls0 ls10 ls20 ls30];
teta0=interp1(ft2,ft1,lw2);
Mv=Pvtr*zp*Av*0.001;
kr1=[-30 -20 -10 0 10 20 30 40 50 60 70 80 90];
N1=[-ls30 -ls20 -ls10 ls0 ls10 ls20 ls30 ls40 ls50 ls60 ls70 ls80 ls90];
M3=[lw2 lw2 lw2 lw2 lw2 lw2 lw2 lw2 lw2 lw2 lw2 lw2 lw2];
% plot(kr1,N1,'x');
Fp1 = griddedInterpolant(kr1,N1);
funp1 = @(t) Fp1(t);
% plot(kr1, funp1(kr1));
bp1 = integral(funp1, teta0, 50);
% plot(kr1,M3,'x');
Fp2 = griddedInterpolant(kr1,M3);
funp2 = @(t) Fp2(t);
% plot(kr1, funp2(kr1));
bp2 = integral(funp2, teta0, 50);
b=bp1-bp2;
% tetar=-24;
kam=0.98;
G1=[2.4 2.5 2.6 2.7 2.8 2.9 3.0 3.1 3.2 3.3 3.4 3.5];
G2=[1 0.98 0.96 0.95 0.93 0.91 0.90 0.88 0.86 0.84 0.82 0.8];
EX1=interp1(G1,G2,x(2)/x(3));
G3=[0.45 0.5 0.55 0.60 0.65 0.7 ];
G4=[0.75 0.82 0.89 0.95 0.97 1];
EX2=interp1(G3,G4,x(5));
r=0.73+0.6*(kz*x(4)-x(3))/x(3);
ctk=0.373+0.023*x(2)/x(3)-0.043*x(1)/100;
Tk=2*ctk*x(2)/sqrt(h0);
G5=[4 5 6 7 8 10 12 14 16 18 20];
G6=[0.1 0.1 0.1 0.098 0.093 0.079 0.065 0.053 0.044 0.038 0.035];
S=interp1(G5,G6,Tk);
tetar=109*kam*EX1*EX2*sqrt(r*S);
ap1=integral(funp1, tetar, 0);
ap2=integral(funp2, tetar, teta0);
ap3=integral(funp1, 0, teta0);
a=abs(ap1)+ap2-ap3;
k=b/a;
%%% Imitation
dem=0;
for kk=1:100
dem=dem+1;
D=D1;
C=404;
M_R=normrnd(300,50,[1,15])
R1=M_R(randperm(numel(M_R),1))
zt=Pt;
P=x(8);
Q=0;
Q2=0;
% ??????
nmonth1 = 31;
month1 = 1.15*rand(nmonth1,1)+0.09;
hw_day1 = randi([1,31]);
month1(hw_day1) = 2.25*rand(1)+1.25;
hw_1=month1.';
% ????????
nmonth2=28;
month2=1.15*rand(nmonth2,1)+0.09;
hw_2=month2.';
% ????
nmonth3=31;
month3=1.15*rand(nmonth3,1)+0.09;
hw_3=month3.';
% ??????
nmonth4=30;
month4 = 1.15*rand(nmonth4,1)+0.09;
hw_day4 = randi([1,31]);
month4(hw_day4) = 2.25*rand(1)+1.25;
hw_4=month4.';
% ???
nmonth5=31;
month5 = 1.15*rand(nmonth5,1)+0.09;
daystorm5=2;
hw_day5 = zeros(2,1);
while hw_day5(1) == hw_day5(2);
for ii=1:daystorm5
hw_day5(ii) = randi([1,nmonth5]);
end
clear ii
end
for ii=1:daystorm5
month5(hw_day5(ii)) = 2.25*rand(1)+1.25;
end
hw_5=month5.';
%????
nmonth6=30;
month6 = 1.15*rand(nmonth6,1)+0.09;
daystorm6=2;
hw_day6 = zeros(2,1);
while hw_day6(1) == hw_day6(2);
for ii=1:daystorm6
hw_day6(ii) = randi([1,nmonth6]);
end
clear ii
end
for ii=1:daystorm6
month6(hw_day6(ii)) = 2.25*rand(1)+1.25;
end
hw_6=month6.';
%????
nmonth7=31;
month7 = 1.15*rand(nmonth7,1)+0.09;
daystorm7=3;
hw_day7 = zeros(3,1);
while hw_day7(1) == hw_day7(2);
for ii=1:daystorm7
hw_day7(ii) = randi([1,nmonth7]);
end
clear ii
end
while hw_day7(1) == hw_day7(3);
for ii=1:daystorm7
hw_day7(ii) = randi([1,nmonth7]);
end
clear ii
end
while hw_day7(2) == hw_day7(3);
for ii=1:daystorm7
hw_day7(ii) = randi([1,nmonth7]);
end
clear ii
end
for ii=1:daystorm7
month7(hw_day7(ii)) = 2.25*rand(1)+1.25;
end
hw_7=month7.';
%??????
nmonth8=31;
month8 = 1.15*rand(nmonth8,1)+0.09;
daystorm8=3;
hw_day8 = zeros(3,1);
while hw_day8(1) == hw_day8(2);
for ii=1:daystorm8
hw_day8(ii) = randi([1,nmonth8]);
end
clear ii
end
while hw_day8(1) == hw_day8(3);
for ii=1:daystorm8
hw_day8(ii) = randi([1,nmonth8]);
end
clear ii
end
while hw_day8(2) == hw_day8(3);
for ii=1:daystorm8
hw_day8(ii) = randi([1,nmonth8]);
end
clear ii
end
for ii=1:daystorm8
month8(hw_day8(ii)) = 2.25*rand(1)+1.25;
end
hw_8=month8.';
%????????
nmonth9=30;
month9 = 1.15*rand(nmonth9,1)+0.09;
daystorm9=3;
hw_day9 = zeros(3,1);
while hw_day9(1) == hw_day9(2);
for ii=1:daystorm9
hw_day9(ii) = randi([1,nmonth9]);
end
clear ii
end
while hw_day9(1) == hw_day9(3);
for ii=1:daystorm9;
hw_day9(ii) = randi([1,nmonth9]);
end
clear ii
end
while hw_day9(2) == hw_day9(3);
for ii=1:daystorm9
hw_day9(ii) = randi([1,nmonth9]);
end
clear ii
end
for ii=1:daystorm9
month9(hw_day9(ii)) = 2.25*rand(1)+1.25;
end
hw_9=month9.';
%???????
nmonth10=31;
month10 = 1.15*rand(nmonth10,1)+0.09;
daystorm10=3;
hw_day10 = zeros(3,1);
while hw_day10(1) == hw_day10(2);
for ii=1:daystorm10
hw_day10(ii) = randi([1,nmonth10]);
end
clear ii
end
while hw_day10(1) == hw_day10(3);
for ii=1:daystorm10
hw_day10(ii) = randi([1,nmonth10]);
end
clear ii
end
while hw_day10(2) == hw_day10(3);
for ii=1:daystorm10
hw_day10(ii) = randi([1,nmonth10]);
end
clear ii
end
for ii=1:daystorm10
month10(hw_day10(ii)) = 2.25*rand(1)+1.25;
end
hw_10=month10.';
% ??????
nmonth11 = 30;
month11 = 1.15*rand(nmonth11,1)+0.09;
hw_day11 = randi([1,31]);
month11(hw_day11) = 2.25*rand(1)+1.25;
hw_11=month11.';
% ???????
nmonth12 = 31;
month12 = 1.15*rand(nmonth12,1)+0.09;
hw_day12 = randi([1,31]);
month12(hw_day12) = 2.25*rand(1)+1.25;
hw_12=month12.';
hw=[hw_1 hw_2 hw_3 hw_4 hw_5 hw_6 hw_7 hw_8 hw_9 hw_10 hw_11 hw_12];
% ???????? ???? ?? ???????
catch1= [3.9 3.6 3.8 3.7 3.8 2.99 3.8 3.9 3.965 4.2 4.26 2.9 3.85 3.7 4.41 3.46 3.58 4.57 3.2 3.59 3.42 3.84 3.96 3.57 3.54 4.43 3.28 4.36 3.85 4.28 3.95];
b_1=mean(catch1);
s_1=std(catch1);
q_1 = normrnd(4.0,s_1,[1,31]);
catch2=[3.92 3.98 3.69 3.64 3.57 4.86 3.89 4.67 3.85 4.12 4.13 3.69 3.92 3.64 3.93 4.58 4.38 4.56 4.65 4.75 3.84 4.28 4.37 4.65 4.28 4.54 4.69 4.38];
b_2=mean(catch2);
s_2=std(catch2);
q_2 = normrnd(4.1,s_2,[1,28]);
catch3=[5.46 4.77 4.82 5.58 4.75 3.98 4.66 4.88 3.95 4.97 4.62 4.58 4.76 4.89 4.86 4.95 5.35 5.52 4.68 4.89 5.95 3.96 4.65 4.36 4.86 3.94 4.84 4.46 4.18 4.76 4.68];
b_3=mean(catch3);
s_3=std(catch3);
q_3 = normrnd(4.2,s_3,[1,31]);
catch4=[4.59 5.10 5.23 5.46 5.28 4.95 4.68 5.36 5.48 5.65 4.28 5.68 5.37 5.85 5.95 4.12 4.58 4.31 4.68 5.36 5.41 4.58 4.68 4.95 5.36 5.75 5.63 5.39 5.48 5.36];
b_4=mean(catch4);
s_4=std(catch4);
q_4 = normrnd(5.1,s_4,[1,30]);
catch5=[5.15 4.85 4.95 5.65 5.89 6.28 6.35 4.95 6.02 6.35 5.89 5.78 5.63 5.94 6.35 6.18 6.75 7.05 5.98 5.04 6.15 6.38 6.54 5.87 7.36 5.84 7.25 6.35 6.25 5.98 5.46];
b_5=mean(catch5);
s_5=std(catch5);
q_5 = normrnd(5.8,s_5,[1,31]);
catch6=[5.84 5.76 5.68 5.38 5.25 5.19 6.25 6.34 7.42 7.57 5.98 6.84 7.05 7.38 6.86 5.89 5.43 5.95 5.96 6.58 5.13 5.38 5.69 7.84 7.26 7.36 7.21 6.89 6.36 7.16];
b_6=mean(catch6);
s_6=std(catch6);
q_6 = normrnd(6.0,s_6,[1,30]);
catch7=[6.98 7.35 6.84 6.95 6.86 6.35 5.45 5.68 5.36 5.98 5.87 6.58 7.35 7.48 7.31 7.29 6.35 7.75 6.98 6.38 6.19 5.64 5.84 7.05 6.89 5.68 7.36 6.35 6.48 6.36 5.36];
b_7=mean(catch7);
s_7=std(catch7);
q_7 = normrnd(6.4,s_7,[1,31]);
catch8=[4.86 5.06 5.34 4.68 4.38 3.89 3.68 3.78 5.42 4.39 4.68 3.57 5.69 4.36 5.39 4.26 4.16 3.95 3.68 3.46 5.39 5.48 5.68 5.31 4.58 5.36 5.36 5.34 5.16 5.27 5.13];
b_8=mean(catch8);
s_8=std(catch8);
q_8 = normrnd(5.4,s_8,[1,31]);
catch9=[4.78 4.16 4.25 4.38 4.69 4.58 4.67 4.23 5.25 5.26 4.98 4.67 4.39 4.68 4.34 4.67 4.85 5.64 4.92 5.39 4.62 5.48 4.13 4.69 4.36 4.30 4.28 3.98 4.71 4.26 ];
b_9=mean(catch9);
s_9=std(catch9);
q_9 = normrnd(5.1,s_9,[1,30]);
catch10=[5.16 5.38 5.24 5.34 5.18 5.39 4.98 4.89 5.32 5.36 5.38 5.29 5.42 5.51 5.25 5.26 5.23 5.34 5.18 5.19 4.86 4.97 4.65 5.01 4.87 4.36 5.21 5.26 5.36 5.08 4.65];
b_10=mean(catch10);
s_10=std(catch10);
q_10 = normrnd(5.4,s_10,[1,31]);
catch11=[5.29 5.68 5.97 6.25 6.35 6.89 6.38 6.85 6.79 6.75 7.24 7.36 7.19 7.28 6.25 6.38 5.98 5.67 6.48 6.35 6.15 6.56 6.68 6.95 7.06 7.25 7.7 6.89 6.38 6.82 ];
b_11=mean(catch11);
s_11=std(catch11);
q_11 = normrnd(6.8,s_11,[1,30]);
catch12=[6.8 6.3 6.4 6.5 5.89 5.98 6.65 6.85 6.87 6.95 6.34 6.85 6.48 7.26 5.39 5.98 6.25 6.28 6.47 5.85 7.26 7.71 7.25 7.65 7.26 6.48 6.68 6.85 7.04 5.89 5.65 ];
b_12=mean(catch12);
s_12=std(catch12);
q_12 = normrnd(6.4,s_12,[1,31]);
q=[q_1 q_2 q_3 q_4 q_5 q_6 q_7 q_8 q_9 q_10 q_11 q_12];
hpr=1.25;
tpod=1;
tneo=2;
t_i=0;
tperi=R1*1.6/(x(9)*1.82)*1.1/24;
tst_i=1;
day=0;
k_i=0;
tu=0;
tcb=tpod+tneo;
j = 1;
z = tperi*qg1*D^(1/3)*x(9)^4.5/C*k1*(1+am);
for day_in_ocean=1:365;
day=day+1;
t_i=day_in_ocean+tu+k_i+tperi+tcb;
if hw(fix(t_i))<hpr;
Q=q(fix(t_i))+Q;
z=z+qg2*(k1)*(1+am)*i1*i2*D^(1/3)*x(9)^4.5/C;
else
tu=tu+tst_i;
z=z+qg2*k1*(1+am)*i1*i2*D^(1/3)*x(9)^4.5/C*0.85;
end
if (z/zt>0.6) || (Q>P);
disp(Q);
j = j+1;
Q_matrix(j)= Q;
z = tperi*qg1*D^(1/3)*x(9)^4.5/C*k1*(1+am);
day=0;
Q2=Q2+Q;
Q=0;
k_i=k_i+2*tperi+tcb;
end
if t_i>285;
% Q1=Q2+Q;
disp(k_i);
break;
end
end
A=Q_matrix;
text = 'Sum of catching fish of each trip';
disp(text);
disp(num2str(Q_matrix));
text = 'Sum of catching fish in a year:';
disp(text);
disp(num2str(Q2));
A(1)=[];
M=mean(A);
Std=std(A);
nr=length(A);
M_txr = normrnd(4400,200,[1,nr]);
dokhod1=0;
for i=1:length(A)
dokhod1=dokhod1+A(i)*M_txr(i);
end
dokhod=dokhod1;
Upop(kk)=dokhod
% x=77;
B=[35 37 39 41 43 45 47 49 50 51 53 55 57 59 61 63 65 67 69 71 73 75 77 79 81 83 85];
for i=1:length(B);
giatri= normpdf(B(i), M, Std);
giatri_matrix(i)=giatri;
end
P=giatri_matrix;
g=sum(P);
pd = makedist('Normal','mu',M,'sigma',Std);
x = 30:85;
y = pdf(pd, x);
end
%%% kinh te
Ks=1.185*10^6;
Tvnhe=12*5+365*0.05;
Tek=365-Tvnhe;
np=Tek/Avt
z1=(250*12+10*Avt+100)*np;
M_txt=normrnd(840,15,[1,15]);
txt=M_txt(randperm(numel(M_txt),1));
txm=2860;
z2=0.9*txt*Pt+0.1*txm*Pt;
z3=x(8)*15;
tziv=750;
z4=0.02*x(8)*tziv;
z5=12*Ks/100;
z6=z5*0.05;
z7=10*Avt*np;
zobsi=z1+z2+z3+z4+z5+z6+z7;
Pr_i=Upop-zobsi*15
E_i=Pr_i/Ks
C_i=1./E_i
M_Ci=mean(C_i);
Std_Ci=std(C_i);
% Q_ii=[1 2 3 4 5 6 7 8 9 10];
Q_ii=[1 1.5 2 2.5 3 3.5 4 4.5 5 5.5 6 6.5 7];
figure(1) %%% this FIGURE
for i=1:length(Q_ii)
Probability= normpdf(Q_ii(i), M_Ci, Std_Ci);
Probability_matrix(i)=Probability;
end
Prob=Probability_matrix(i);
giatri2=4.83;
Pro_giatri2= normpdf(giatri2, M_Ci, Std_Ci);
pd_Ci = makedist('Normal','mu',M_Ci,'sigma',Std_Ci);
xx = 1:0.1:7;
yy = pdf(pd_Ci, xx);
hsrok=histogram(C_i,'Normalization','pdf');
hsrok.NumBins = 10;
hsrok.BinEdges = [2.5:6.5]
hsrok.BinWidth = 0.25;
hsrok.FaceColor='y';
hold on;
plot(xx,yy, 'linewidth',1);
xlabel('Payback')
ylabel('frequency')
grid on;
% Loi nhuan
figure(2) %%% this FIGURE
Pr_i=Pr_i/1000;
M_Pri=mean(Pr_i);
Std_Pri=std(Pr_i);
giatri1=245.417
Prob_giatri1= normpdf(giatri1, M_Pri, Std_Pri)*10;
Pr_ii=[165000 180000 195000 210000 225000 240000 255000 270000 285000 300000 315000 330000 345000 360000 375000 390000]/1000;
for i=1:length(Pr_ii)
Prob_Pr= normpdf(Pr_ii(i), M_Pri, Std_Pri);
Prob_Pr_matrix(i)=Prob_Pr;
end
Prob_Pr1=Prob_Pr_matrix;
pd_Pri = makedist('Normal','mu',M_Pri,'sigma',Std_Pri);
xxx = 100:0.1:450;
yyy = pdf(pd_Pri, xxx);
hpri=histogram(Pr_i,'Normalization','pdf');
hpri.NumBins = 12;
hpri.FaceColor='y';
hold on
plot(xxx,yyy,'linewidth',1);
xlabel('Profit')
ylabel('frequency')
grid on;
xacsuat1=Pro_giatri2;
xacsuat1=num2str(xacsuat1);
sr_11=min(C_i) % thoi gian thu hoi von
f=sr_11 % gia tri ham muc tieu
end
function [c, ceq] = nonlcond(x)
global Wgr h0 ls k Tk D1 sr_11
c=[90-Wgr
0.35-h0
0.25-ls
1-k
4.6-Tk
3-sr_11];
ceq = [1.002*x(1)*x(2)*x(3)*x(5)*1.025-D1];
end
disp(history.x);
disp(history.fval);
disp(searchdir)
disp(figure(2))
disp(C_i); % This C_i=[]
pkll=C_i
end

请先登录，再回答此问题。

回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Solver Outputs and Iterative Display 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by