edit matlab java code
Show older comments
I1max = 3; %最大部品1の在庫量(3)
I1min = 0; %最小部品1の在庫量(0)
I2max = 3; %最大部品2の在庫量(3)
I2min= 0; %最小部品2の在庫量(0）
I3max =3; %最大配送中の在庫量(3)
I3min= 0; %最小配送中の在庫量(0)
Jmax = 6; %最大製品の在庫量(6)
Jmin = -1; %最小製品の在庫量(-1)
C = 5; %生産(工程)能力(5)
c1 = 4; %部品1の再生産費
c2 = 4; %部品2の再生産費
u2u=0;%回収業者2UPの状態
u2d=1;%回収業者2DOWNの状態
f=1;%回収業者2ＵＰからUPの確率;
g=1;%回収業者2ＤＯＷＮからUPの確率；
%需要の決定
Q=2; %分散
RD=2; %平均需要
Dnmin=RD-(Q/2);
Dnmax=RD-(Q/2)+Q;
N1=nchoosek(Q,0);
N2=nchoosek(Q,1);
N3=nchoosek(Q,2);
D(1)=N1*(1/2)^Q;
D(2)=N2*(1/2)^Q;
D(3)=N3*(1/2)^Q;
%決定kの最大値を算出
Dmlgh = (I1max+I1min+1)*(I3max+I2max-1)*(I2max-I2min-1)*(Jmax+(-Jmin)+1)*(u2u+u2d+1);
L=Dmlgh;
for i1=I1min:I1max
fori2=I2min:I2max
fori3=I3min:I3max-i2
Y=myfunction(i2,i3);
forj=Jmin:Jmax
for u2=u2u:u2d
l=(u2-u2u+1)+(j-Jmin)*(u2d-u2u+1)+(Jmax-Jmin+1)*Y*(u2d-u2u+1)+(Jmax-Jmin+1)*(I3max+I2max-1)*(I2max-I2min-1)*(i1-I1min)*(u2d-u2u+1);
k1max(l)=I1max-i1; %k1決定空間
if u2==0
k2max(l)=(I2max-(i2+i3)); %k2決定空間
else
k2max(l)=0;%k2決定空間
end
if(Jmax-j+Dnmin)>=i2+i1
ifi2+i1>=C %Cは工程（生産）能力
k3max(l)=C;
else
k3max(l)=i2+i1;
end
end
if(Jmax-j+Dnmin)<=i2+i1
ifJmax-j+Dnmin>=C
k3max(l)=C;
else
k3max(l)=Jmax-j+Dnmin;
end
end
ifk1max(l)<0
k1max(l)=0;
end
ifk1max(l)>I1max
k1max(l)=I1max;
end
ifk2max(l)>I2max
k2max(l)=I2max;
end
ifk2max(l)<0
k2max(l)=0;
end
if k3max(l)<0
k3max(l)=0;
end
ifJmax-j+Dnmin>i2+i1
ifk3max(l)<0
k3max(l)=0;
end
ifk3max(l)>I1max+I2max
k3max(l)=I1max+I2max;
end
end
end
end
end
end
end
decmax=0;
%推移確率pijの計算の準備段階
for i1=I1min:I1max
fori2=I2min:I2max
fori3=I3min:I3max-i2
Y=myfunction(i2,i3);
forj=Jmin:Jmax
for u2=u2u:u2d
l=(u2-u2u+1)+(j-Jmin)*(u2d-u2u+1)+(Jmax-Jmin+1)*Y*(u2d-u2u+1)+(Jmax-Jmin+1)*(I3max+I2max-1)*(I2max-I2min-1)*(i1-I1min)*(u2d-u2u+1); %ｔ期
fork=0:k1max(l)
fork13=0:k2max(l)
fork21=0:k3max(l)
dec=(k21+1)+(k3max(l)+1)*(k13)+(k3max(l)+1)*(k2max(l)+1)*(k);%決定
ifdec>decmax
decmax=dec;
end
%推移確率pijの値がセットされるところに初期値 0 をセット
for ll=1:Dmlgh
P(l,ll,dec)=0;
end
for U2=u2u:u2d
a=u2*2+U2;
if a<1
p2=f;
elseif a==1
p2=(1-f);
elseif a==2;
p2=g;
else
p2=(1-g);
end
if a==0
A=1;
else
A=0;
end
ford=Dnmin:Dnmax
if k21<=i2
I1=round(i1+k);
I3=round((k13)*A);
I2=round(i2+i3*A-k21);
J=round(j-d+k21);
U=round(U2);
else
I1=round(i1+k-k21+i2);
I3=round((k13)*A);
I2=round(i3*A);
J=round(j-d+k21);
U=round(U2);
end
if J>Jmax
J=Jmax; end
if J<Jmin
J=Jmin; end
ifI1>I1max
I1=I1max; end
ifI1<I1min
I1=I1min; end
ifI2>I2max
I2=I2max; end
ifI2<I2min
I2=I2min; end
ifI3>I3max
I3=I3max; end
ifI3<I3min
I3=I3min; end
T=myfunction(I2,I3);
ll=(U-u2u+1)+(J-Jmin)*(u2d-u2u+1)+(Jmax-Jmin+1)*T*(u2d-u2u+1)+(Jmax-Jmin+1)*(I3max+I2max-1)*(I2max-I2min-1)*(I1-I1min)*(u2d-u2u+1);
P(l,ll,dec)=P(l,ll,dec)+D(d)*p2;%推移確率
end
end
end
end
end
end
end
end
end
end
a1=5; % 購入費(3)
a2=1; % 購入費
hi1=1; % 部品保管費用(1)
hi2=1; % 部品保管費用(1)
hj=2; % 製品保管費用(2)
cb=20; % 機会損失費(20)
for l=1:Dmlgh
fordec=1:decmax
R(l,dec)=10000;
end
end
fori1=I1min:I1max
fori2=I2min:I2max
fori3=I3min:I3max-i2
Y=myfunction(i2,i3);
forj=Jmin:Jmax
for u2=u2u:u2d
l=(u2-u2u+1)+(j-Jmin)*(u2d-u2u+1)+(Jmax-Jmin+1)*Y*(u2d-u2u+1)+(Jmax-Jmin+1)*(I3max+I2max-1)*(I2max-I2min-1)*(i1-I1min)*(u2d-u2u+1);
fork=0:k1max(l)
fork13=0:k2max(l)
for k21=0:k3max(l)
dec=(k21+1)+(k3max(l)+1)*(k13)+(k3max(l)+1)*(k2max(l)+1)*(k);
ifdec>decmax
decmax=dec;
end
ford=Dnmin:Dnmax
if k21<=i2
r=a1*k+a2*k13*(f^2)+c2*k21+hi1*i1+hi2*i2+hj*max(0,j)+cb*max(0,d-j);
else
r=a1*k+a2*k13*(f^2)+c1*(k21-i2)+c2*i2+hi1*i1+hi2*i2+hj*max(0,j)+cb*max(0,d-j);
end
R(l,dec)=r;
end
end
end
end
end
end
end
end
end
P(l,ll,dec)
R(l,dec)
[policy,average_reward,cpu_time]=mdp_relative_value_iteration(P,R);%MDPtoolboxを使う
%最適政策
Q=zeros(L,L);
for t=1:L
opk=policy(t);
k1m=k1max(t);
k2m=k2max(t);
k3m=k3max(t);
Q(t,:)=P(t,:,opk);
[part1(t),part2(t),product(t)]=E_policy(opk,k1m,k2m,k3m);
End
function[par1,par2,pro] = E_policy(Po,k1m,k2m,k3m)
for par1=0:k1m
for par2=0:k2m
for pro=0:k3m
if(pro)+(k3m+1)*(par2)+(k3m+1)*(k2m+1)*(par1)==Po
disp(par1)
disp(par2)
disp(pro)
return;
end
end
end
end
function y=myfunction(g,f)
if g==2
y=7+f;
elseif g==3
y=9+f;
else
y=4*g+f;
end
1 Comment
Answers (0)
Categories
Find more on 起動と終了 in Help Center and File Exchange
on 8 Oct 2019
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
