Borrar filtros
Borrar filtros

Imaginary Component in Total Profit Calculation

7 visualizaciones (últimos 30 días)
Shaf
Shaf el 2 de Sept. de 2024
Editada: John D'Errico el 2 de Sept. de 2024
Hi everyone,
I'm working on a three-echelon inventory model involving a farmer, processor, and retailer, and I'm encountering an issue while calculating the total profit. The result I get is 2.6631e+04 - 1.0581e+05i, which includes an unexpected imaginary component.
The model has three decision variables:
1. n (shipment frequency),
2. g (preservation technology cost),
3. x (replenishment cycle of the retailer (days)).
Here’s a brief overview of the algorithm I'm using:
1. Initialize n=1
2. Initialize parameters x(1) = 1 and g(1) = 0.0007.
3. Calculate x(2) using x(1) and g(1).
4. Calculate g(2) using x(2) and g(1).
5. Repeat the calculations until changes in x and g between iterations are minimal.
6. Compute total profit.
7. Set n=2
8. Repeat step 2-7
9. Compare the total profit for n=1 and n=2, if the total profit increases, continue the algorithm.
Here’s the MATLAB code I’m using:
a = 80;%Scale parameter of the demand rate
b = 20;%Shape parameter of demand rate
bk = 0.2;%mortality rate
c = 0.2;%carbon tax
cf = 0.1;%feeding cost
cg = 0.01; %carbon emission
f = 0.9;%survival rate
hp = 0.2; %holding cost processor/manufacture
hr = 0.1; %holding cost retailer
hs = 0.07; %carbon emission
ht = 0.05; %carbon emission
Ik = 0.1; %interest charge (credit payment)
j = 0.001; %carbon emission
k = 6.87; %Asymptotic weight of the items
Kf = 500; %ordering cost farmer
Kj = 0.3; %carbon emission
Km = 0.5; %carbon emission
Kp = 100; %setup cost manufacture
Kr = 100; %ordering cost retailer (100-1000)
Ks = 1; %carbon emission
l = 0.3; %interest earned (credit payment)
L = 7; %uproduct lifetime/lifespan (day)
m = 0.2; %mortaility cost
M = 15; %credit payment period (day)
ml = 0.01; %carbon emission
pf = 0.7; %purchasing cost farmer
pp = 2; %purchasing cost processor/manufacture
pr = 6; %purchjasing cost retailer
pt = 0.2; %carbon emission
Pv = 1;%procurement cost farmer
pw = 0.5; %carbon emission
q = 0.1; %Exponential growth rate of the items (0.1-1)
R = 250; %production rate
S = 100; %transportation cost retailer
Sr = 1; %carbon emission
Sv = 100; %transportation cost processor/manufacture
T = 0.5; %carbon emission
v = 0.15; %production cost
vp = 100; %carbon emission
w = 0.064; %initial weight
pq = 0.5; %carbon emission
ww = 5; %weight target (1-100)
y = 335; %Scale factor for the preservation technology function (100-1000, 335)
z = 10; %constant of integration (10-100)
Tf=-log(((k/ww)-1)/z)/q; %T farmer
Stop = 0
n=1;
fy=2;
JTP(fy-1)=0;
while Stop==0
x(1)= 0.6;
g(1)= 0.0007;
xa = 2*(Kr + c*Ks) + 2*(Kf + c*(Kj + Km) + Kp)/n + 2^(b/(-1+b))*c*(pf + pp)*(-(a*(-1+b)*g(1)*x(1)^2*y)/(L + g(1)*L*y))^(1/(1-b))+(2^((-1+2*b)/(-1+b))*c*(pf+pp)*(-(a*(-1+b)*g(1)*x(1)^2*y)/(L + g(1)*L*y))^(1/(1-b)))/(-1+b);
xb = 2^(b/(-1+b))*pr*(-(a*(-1+b)*g(1)*x(1)^2*y)/(L + g(1)*L*y))^(1/(1-b))-(2^((-1+2*b)/(-1+b))*pr*(-(a*(-1+b)*g(1)*x(1)^2*y)/(L + g(1)*L*y))^(1/(1-b)))/(-1+b)+(2^(b/(-1+b))*(c*v + vp)*(-(a*(-1+b)*g(1)*x(1)^2*y)/(L + g(1)*L*y))^(1/(1-b)))/n;
xc = (2^((-1+2*b)/(-1+b))*(c*v + vp)*(-(a*(-1+b)*g(1)*x(1)^2*y)/(L + g(1)*L*y))^(1/(1-b)))/((-1+b)*n)+(2^(b/(-1+b))*(Pv + c*pw)*w*(-(a*(-1+b)*g(1)*x(1)^2*y)/(L + g(1)*L*y))^(1/(1-b)))/(f*ww);
xd = (2^((-1+2*b)/(-1+b))*(Pv+c*pw)*w*(-(a*(-1+b)*g(1)*x(1)^2*y)/(L + g(1)*L*y))^(1/(1-b)))/((-1+b)*f*ww)+4^(1/(-1+b))*(g(1) + hp + c*(ht + j))*(-1+n)*(-(a*(-1+b)*g(1)*x(1)^2*y)/(L + g(1)*L*y))^(-2/(-1+b));
xe = (2^(2*b/(-1+b))*(g(1) + hp + c*(ht + j))*(-1+n)*(-(a*(-1+b)*g(1)*x(1)^2*y)/(L + g(1)*L*y))^(-2/(-1+b)))/(-1+b)+(4^(1/(-1+b))*(g(1)+hp+c*(ht+j))*(-(a*(-1+b)*g(1)*x(1)^2*y)/(L+g(1)*L*y))^(-2/(-1+b)))/R;
xf = (2^(2*b/(-1+b))*(g(1) + hp + c*(ht + j))*(-(a*(-1+b)*g(1)*x(1)^2*y)/(L + g(1)*L*y))^(-2/(-1+b)))/((-1+b)*R)-2*l*pr*(2^(1/(-1+b))*(-M+x(1))*(-(a*(-1+b)*g(1)*x(1)^2*y)/(L + g(1)*L*y))^(1/(1-b)))+(2^(1/(-1+b))*l*pr*4*L*M*(-(a*(-1+b)*g(1)*x(1)^2*y)/(L + g(1)*L*y))^(1/(1-b)))/((-1+b)*L*(1 + g(1)*y));
xg = (2^(b/(-1+b)) * k*(cf*f + c*cg*f + bk*m + bk*c*ml) * (-(a*(-1+b)*g(1)*x(1)^2*y)/(L + g(1)*L*y))^(1/(1-b))*(q*Tf-log(1+z)+log(1+exp(-q*Tf)*z)))/(f*q*ww)+(2^((-1+2*b)/(-1+b))*k*(cf*f+c*cg*f+bk*m+bk*c*ml)*x(1)*(-(a*(-1+b)*g(1)*x(1)^2*y)/(L + g(1)*L*y))^(1/(1-b))*(q*Tf-log(1+z)+log(1+exp(-q*Tf)*z)))/((-1+b)*f*q*ww);
xh = (1/(2*x(1)^3))*x(1)*(xa-xb+xc+xd+xe+xf+xg);
xj = -(((8*(g(1)+hr+c*(hs+j))*(-(a*(-1+b)*g(1)*x(1)^2*y)/(L+g(1)*L*y))^(1/(1-b)))/(-1+b)+(2^(1/(-1+b))*l*pr*(2*L*(2*M+(-3+b))*(-(a*(-1+b)*g(1)*x(1)^2*y)/(L+g(1)*L*y))^(1/(1-b)) + 2*(-3+b)*L*(-(a*(-1+b)*g(1)*x(1)^2*y)/(L+g(1)*L*y))^(1/(1-b))))/((-1+b)*L*(1+g(1)*y))+x(1)*(2^(-2+b/(-1+b))*g(1)*y*((a*(1-b)*g(1)*x(1)^2*y)/(L+g(1)*L*y))^(b/(1-b)))/(L+g(1)*L*y)+x(1)^2*(2^(1/(-1+b))*l*pr*g(1)*y*(-3+b)*((a*(1-b)*g(1)*x(1)^2*y)/(L+g(1)*L*y))^(b/(1-b)))/((-1+b)*L*(1+g(1)*y)))/(2*x(1)));
xk = (4*(c*Sr+Sv)/n^2+4*(S+c*T));
x(2) = ((xk+xh)/xj)^(1/2);
ga(2)=2^(1/(-1+b))*L*n*(Pv+c*pw)*q*R*w*(1+g(1)*y)*(-(a*(-1+b)*x(2)^2*y)/(L+g(1)*L*y))^(1/(1-b));
gb(2)=2^(1/(-1+b))*cf*L*n*(pf+pp)*q*R*ww*(1+g(1)*y)*(-(a*(-1+b)*x(2)^2*y)/(L+g(1)*L*y))^(1/(1-b));
gc(2)=-2^(1/(-1+b))*f*L*n*pr*q*R*ww*(1+g(1)*y)*(-(a*(-1+b)*x(2)^2*y)/(L+g(1)*L*y))^(1/(1-b))+2^(1/(-1+b))*f*L*q*R*(c*v+vp)*ww*(1+g(1)*y)*(-(a*(-1+b)*x(2)^2*y)/(L+g(1)*L*y))^(1/(1-b));
gd(2)=2*f*(g(1)+hr+c*(hs+j))*L*n*q*R*ww*x(2)*(1+g(1)*y)*(-(a*(-1+b)*x(2)^2*y)/(L+g(1)*L*y))^(1/(1-b))-2*(-1+b)*f*g(1)*L*n*q*R*ww*x(2)*(1+g(1)*y)^2*(-(a*(-1+b)*x(2)^2*y)/(L+g(1)*L*y))^(1/(1-b));
ge(2)=2^((3-2*b)/(-1+b))*f*l*n*pr*q*R*ww*(4*L*(M-x(2))*(-(a*(-1+b)*x(2)^2*y)/(L+g(1)*L*y))^(1/(1-b))-2*g(1)*y*(x(2)^3*((a*(1-b)*x(2)^2*y)/(L+g(1)*L*y))^(b/(1-b)) - g(1)^(1/(1-b))*2*L*M*(-(a*(-1+b)*x(2)^2*y)/(L+g(1)*L*y))^(1/(1-b))+2*L*x(2)*(-(a*(-1+b)*x(2)^2*y)/(L+g(1)*L*y))^(1/(1-b))))+2^(1/(-1+b))*k*L*(cf*f+c*cg*f+bk*m+bk*c*ml)*n*R*(1+g(1)*y)*(-(a*(-1+b)*x(2)^2*y)/(L+g(1)*L*y))^(1/(1-b))*(q*Tf-log(1+z) + log(1+exp(-q*Tf)*z));
gi(2)=-4^(1/(-1+b))*f*(g(1)+hp+c*(ht+j))*L*n*ww*(1+g(1)*y)*(-(a*(-1+b)*x(2)^2*y)/(L+g(1)*L*y))^(-2/(-1+b)) ...
+ 4^(1/(-1+b))*f*(g(1)+hp+c*(ht+j))*L*(1-n)*n*q*R*ww*(1+g(1)*y)*(-(a*(-1+b)*x(2)^2*y)/(L+g(1)*L*y))^(-2/(-1+b)) ...
+ 2^((3-b)/(-1+b))*(-1+b)*f*L*n*ww*(1+g(1)*y)^2*g(1)*(-(a*(-1+b)*x(2)^2*y)/(L+g(1)*L*y))^(-2/(-1+b)) ...
+ g(1)^2*((3-b)/(-1+b))*(-1+b)*f*L*(1-n)*n*q*R*ww*(1+g(1)*y)^2*(-(a*(-1+b)*x(2)^2*y)/(L+g(1)*L*y))^(-2/(-1+b));
g(2)=((ga(2)+gb(2)+gc(2)+gd(2)+ge(2))/gi(2))^(1-b);
i=2;
while abs(x(i)-x(i-1))>0.001 & abs(g(i)-g(i-1))>0.001
i=i+1;
xa(i)=2*(Kr + c*Ks) + 2*(Kf + c*(Kj + Km) + Kp)/n + 2^(b/(-1+b)) * c*(pf + pp) * (-(a*(-1+b)*g(i-1)*x(i-1)^2*y)/(L + g(i-1)*L*y))^(1/(1-b))+(2^((-1+2*b)/(-1+b))*c*(pf+pp)*(-(a*(-1+b)*g(i-1)*x(i-1)^2*y)/(L + g(i-1)*L*y))^(1/(1-b)))/(-1+b);
xb(i)=2^(b/(-1+b))*pr*(-(a*(-1+b)*g(i-1)*x(i-1)^2*y)/(L + g(i-1)*L*y))^(1/(1-b))-(2^((-1+2*b)/(-1+b))*pr*(-(a*(-1+b)*g(i-1)*x(i-1)^2*y)/(L + g(i-1)*L*y))^(1/(1-b)))/(-1+b)+(2^(b/(-1+b))*(c*v + vp)*(-(a*(-1+b)*g(i-1)*x(i-1)^2*y)/(L + g(i-1)*L*y))^(1/(1-b)))/n;
xc(i)=(2^((-1+2*b)/(-1+b))*(c*v + vp)*(-(a*(-1+b)*g(i-1)*x(i-1)^2*y)/(L + g(i-1)*L*y))^(1/(1-b)))/((-1+b)*n)+(2^(b/(-1+b))*(Pv + c*pw)*w*(-(a*(-1+b)*g(i-1)*x(i-1)^2*y)/(L + g(i-1)*L*y))^(1/(1-b)))/(f*ww);
xd(i)=(2^((-1+2*b)/(-1+b))*(Pv+c*pw)*w*(-(a*(-1+b)*g(i-1)*x(i-1)^2*y)/(L + g(i-1)*L*y))^(1/(1-b)))/((-1+b)*f*ww)+4^(1/(-1+b))*(g(i-1) + hp + c*(ht + j))*(-1+n)*(-(a*(-1+b)*g(i-1)*x(i-1)^2*y)/(L + g(i-1)*L*y))^(-2/(-1+b));
xe(i)=(2^(2*b/(-1+b))*(g(i-1) + hp + c*(ht + j))*(-1+n)*(-(a*(-1+b)*g(i-1)*x(i-1)^2*y)/(L + g(i-1)*L*y))^(-2/(-1+b)))/(-1+b)+(4^(1/(-1+b))*(g(i-1)+hp+c*(ht+j))*(-(a*(-1+b)*g(i-1)*x(i-1)^2*y)/(L+g(i-1)*L*y))^(-2/(-1+b)))/R;
xf(i)=(2^(2*b/(-1+b))*(g(i-1) + hp + c*(ht + j))*(-(a*(-1+b)*g(i-1)*x(i-1)^2*y)/(L + g(i-1)*L*y))^(-2/(-1+b)))/((-1+b)*R)-2*l*pr*(2^(1/(-1+b))*(-M+x(i-1))*(-(a*(-1+b)*g(i-1)*x(i-1)^2*y)/(L + g(i-1)*L*y))^(1/(1-b)))+(2^(1/(-1+b))*l*pr*4*L*M*(-(a*(-1+b)*g(i-1)*x(i-1)^2*y)/(L + g(i-1)*L*y))^(1/(1-b)))/((-1+b)*L*(1 + g(i-1)*y));
xg(i)=(2^(b/(-1+b)) * k*(cf*f + c*cg*f + bk*m + bk*c*ml) * (-(a*(-1+b)*g(i-1)*x(i-1)^2*y)/(L + g(i-1)*L*y))^(1/(1-b))*(q*Tf-log(1+z)+log(1+exp(-q*Tf)*z)))/(f*q*ww)+(2^((-1+2*b)/(-1+b))*k*(cf*f+c*cg*f+bk*m+bk*c*ml)*x(i-1)*(-(a*(-1+b)*g(i-1)*x(i-1)^2*y)/(L + g(i-1)*L*y))^(1/(1-b))*(q*Tf-log(1+z)+log(1+exp(-q*Tf)*z)))/((-1+b)*f*q*ww);
xh(i)=(1/(2*x(i-1)^3))*x(i-1)*(xa(i)-xb(i)+xc(i)+xd(i)+xe(i)+xf(i)+xg(i-1));
xj(i)=((8*(g(i-1)+hr+c*(hs+j))*(-(a*(-1+b)*g(i-1)*x(i-1)^2*y)/(L+g(i-1)*L*y))^(1/(1-b)))/(-1+b)+(2^(1/(-1+b))*l*pr*(2*L*(2*M+(-3+b))*(-(a*(-1+b)*g(i-1)*x(i-1)^2*y)/(L+g(i-1)*L*y))^(1/(1-b)) + 2*(-3+b)*L*(-(a*(-1+b)*g(i-1)*x(i-1)^2*y)/(L+g(i-1)*L*y))^(1/(1-b))))/((-1+b)*L*(1+g(i-1)*y))+x(i-1)*(2^(-2+b/(-1+b))*g(i-1)*y*((a*(1-b)*g(i-1)*x(i-1)^2*y)/(L+g(i-1)*L*y))^(b/(1-b)))/(L+g(i-1)*L*y)+x(i-1)^2*(2^(1/(-1+b))*l*pr*g(i-1)*y*(-3+b)*((a*(1-b)*g(i-1)*x(i-1)^2*y)/(L+g(i-1)*L*y))^(b/(1-b)))/((-1+b)*L*(1+g(i-1)*y)))/(2*x(i-1));
xk(i)=(4*(c*Sr+Sv)/n^2+4*(S+c*T));
x(i) = sqrt((xk(i)+xh(i))/xj(i));
ga(i)=2^(1/(-1+b))*L*n*(Pv+c*pw)*q*R*w*(1+g(i-1)*y)*(-(a*(-1+b)*x(i)^2*y)/(L+g(i-1)*L*y))^(1/(1-b));
gb(i)=2^(1/(-1+b))*cf*L*n*(pf+pp)*q*R*ww*(1+g(i-1)*y)*(-(a*(-1+b)*x(i)^2*y)/(L+g(i-1)*L*y))^(1/(1-b));
gc(i)=-2^(1/(-1+b))*f*L*n*pr*q*R*ww*(1+g(i-1)*y)*(-(a*(-1+b)*x(i)^2*y)/(L+g(i-1)*L*y))^(1/(1-b))+2^(1/(-1+b))*f*L*q*R*(c*v+vp)*ww*(1+g(i-1)*y)*(-(a*(-1+b)*x(i)^2*y)/(L+g(i-1)*L*y))^(1/(1-b));
gd(i)=2*f*(g(i-1)+hr+c*(hs+j))*L*n*q*R*ww*x(i)*(1+g(i-1)*y)*(-(a*(-1+b)*x(i)^2*y)/(L+g(i-1)*L*y))^(1/(1-b))-2*(-1+b)*f*g(i-1)*L*n*q*R*ww*x(i)*(1+g(i-1)*y)^2*(-(a*(-1+b)*x(i)^2*y)/(L+g(i-1)*L*y))^(1/(1-b));
ge(i)=2^((3-2*b)/(-1+b))*f*l*n*pr*q*R*ww*(4*L*(M-x(i))*(-(a*(-1+b)*x(i)^2*y)/(L+g(i-1)*L*y))^(1/(1-b))-2*g(i-1)*y*(x(i)^3*((a*(1-b)*x(i)^2*y)/(L+g(i-1)*L*y))^(b/(1-b)) - g(i-1)^(1/(1-b))*2*L*M*(-(a*(-1+b)*x(i)^2*y)/(L+g(i-1)*L*y))^(1/(1-b))+2*L*x(i)*(-(a*(-1+b)*x(i)^2*y)/(L+g(i-1)*L*y))^(1/(1-b))))+2^(1/(-1+b))*k*L*(cf*f+c*cg*f+bk*m+bk*c*ml)*n*R*(1+g(i-1)*y)*(-(a*(-1+b)*x(i)^2*y)/(L+g(i-1)*L*y))^(1/(1-b))*(q*Tf-log(1+z) + log(1+exp(-q*Tf)*z));
gi(i)=4^(1/(-1+b))*f*(g(i-1)+hp+c*(ht+j))*L*n*ww*(1+g(i-1)*y)*(-(a*(-1+b)*x(i)^2*y)/(L+g(i-1)*L*y))^(-2/(-1+b)) ...
+ 4^(1/(-1+b))*f*(g(i-1)+hp+c*(ht+j))*L*(1-n)*n*q*R*ww*(1+g(i-1)*y)*(-(a*(-1+b)*x(i)^2*y)/(L+g(i-1)*L*y))^(-2/(-1+b)) ...
+ 2^((3-b)/(-1+b))*(-1+b)*f*L*n*ww*(1+g(i-1)*y)^2*g(i-1)*(-(a*(-1+b)*x(i)^2*y)/(L+g(i-1)*L*y))^(-2/(-1+b)) ...
+ g(i-1)^2*((3-b)/(-1+b))*(-1+b)*f*L*(1-n)*n*q*R*ww*(1+g(i-1)*y)^2*(-(a*(-1+b)*x(i)^2*y)/(L+g(i-1)*L*y))^(-2/(-1+b));
g(i)=(ga(i)/gi(i))^(1-b);
end
x_opt(n)=x(i)
g_opt(n)=g(i)
%profit for supplier, processor/manufacture,retailer
JTR(fy) = (pr/x_opt(n))*((a*(1-b)*x_opt(n)^2)/(2*L)*((y*g_opt(n))/(1+y*g_opt(n))))^(1/(1-b))-(hr + c*hs + g_opt(n) + c*j)/(x_opt(n) * (2*x_opt(n) * (-(a*(-1+b)*g_opt(n)*x_opt(n)^2*y)/(L + g_opt(n)*L*y)))^(1/(1-b)))-(Kr + c*Ks)/x_opt(n)-(pp + c*pt)/x_opt(n) * ((a*(1-b)*x_opt(n)^2)/(2*L) * ((y*g_opt(n))/(1 + y*g_opt(n))))^(1/(1-b))-(S + c*T)/x_opt(n)^2+(pr/x_opt(n)) * l * ((2^(-2 + b/(-1 + b))*x_opt(n)^3*g_opt(n)*y * ((a*(1-b)*x_opt(n)^2*g_opt(n)*y)/(L + g_opt(n)*L*y))^(b/(1-b)))/(L + g_opt(n)*L*y)-((a*(1-b)*x_opt(n)^2)/(2*L)*((y*g_opt(n))/(1 + y*g_opt(n))))^(1/(1-b))*(M - x_opt(n)));
JTM(fy)= (pp/x_opt(n))*((a*(1-b)*x_opt(n)^2)/(2*L)*((y*g_opt(n))/(1+y*g_opt(n))))^(1/(1-b)) - (Kp+c*Km)/(n*x_opt(n)) - (pf+c*pq)/x_opt(n)*((a*(1-b)*x_opt(n)^2)/(2*L)*((y*g_opt(n))/(1+y*g_opt(n))))^(1/(1-b)) - (hp+c*ht+g_opt(n)+c*j)/(2*R*x_opt(n))*((a*(1-b)*x_opt(n)^2)/(2*L)*((y*g_opt(n))/(1+y*g_opt(n))))^(2/(1-b)) - ((hp+c*ht+g_opt(n)+c*j)*(n-1))/(2*x_opt(n))*((a*(1-b)*x_opt(n)^2)/(2*L)*((y*g_opt(n))/(1+y*g_opt(n))))^(2/(1-b)) - (Sv+c*Sr)/(n*x_opt(n))^2 - ((vp+c*v)/(n*x_opt(n)))*((a*(1-b)*x_opt(n)^2)/(2*L)*((y*g_opt(n))/(1+y*g_opt(n))))^(1/(1-b));
JTS(fy) = (pf/x_opt(n))*((a*(1-b)*x_opt(n)^2)/(2*L)*((y*g_opt(n))/(1+y*g_opt(n))))^(1/(1-b)) - ((Pv+c*pw)*w)/(x_opt(n)*f*ww)*((a*(1-b)*x_opt(n)^2)/(2*L)*((y*g_opt(n))/(1+y*g_opt(n))))^(1/(1-b)) - (cf*f+c*cg*f+m*bk+c*ml*bk)/(x_opt(n)*f*ww)*((a*(1-b)*x_opt(n)^2)/(2*L)*((y*g_opt(n))/(1+y*g_opt(n))))^(1/(1-b))*(k*Tf+k/q*(log(1+z*exp(-q*Tf))-log(1+z)));
%total profit for the whole three eschelon supply chain
JTP(fy)= JTS(fy)+JTM(fy)+JTR(fy);
if JTP(fy)>JTP(fy-1)
n=n+1;
fy=fy+1;
else
Stop=1
end
end
n_optimal=n-1
x_optimal=x_opt(n-1)
g_optimal=g_opt(n-1)
TotProfitProcessor= JTM(fy-1)
TotProfitRetailer = JTR(fy-1)
TotProfitFarmer = JTS(fy-1)
TotProfit=JTP(fy-1)
I’ve already searched through this forum for similar issues, but I still don’t fully understand why the imaginary part is appearing. This imaginary component not only shows up in the total profit but also in the profits for the manufacturer, farmer, and retailer, as well as in the values of the decision variables.
This is my first time using MATLAB, so I’m wondering if there’s anything specific I should check to ensure that no imaginary component appears in the results. Any insights on why this might be happening and how I can resolve it would be greatly appreciated!
Thanks in advance for your help!

Respuestas (2)

Shashi Kiran
Shashi Kiran el 2 de Sept. de 2024
Hi @Shaf,
I have reviewed your code and found the following insights:
When calculating xa, the term:
term = (-(a*(-1+b)*g(i-1)*x(i-1)^2*y)/(L + g(i-1)*L*y))^(1/(1-b))+(2^((-1+2*b)/(-1+b))
the above term simplifies to , this is where an imaginary number appears, and similarly in other calculations of xi.
If you can adjust this part of the implementation, the profit should yield a real value.
Hope this helps.
  4 comentarios
Shaf
Shaf el 2 de Sept. de 2024
Alright, so basically the equation of x is obtained from the derivation of total profit equation. Here’s the total profit equation:
Below is the derivation of that equation:
And here’s the equation for x that I obtained from the derivation. Since it’s quite long, I’ve broken it down into several parts:
Please let me know if anything is unclear or if this isn’t what you were referring to.
John D'Errico
John D'Errico el 2 de Sept. de 2024
Editada: John D'Errico el 2 de Sept. de 2024
@Shashi Kiran is exactly on target here. The problem arises from trying to raise a negative number to a fractional power. This will always result in an imaginary part. Even if the power is a simple one, like 1/3, because the principle cube root is not the one you expect. Try it.
(-1)^(1/3)
ans = 0.5000 + 0.8660i
And while we expect to see this happen:
nthroot(-1,3)
ans = -1
it does not. Worse, when the exponent is a totally general real number with a non-trivial fractional part, you must always get an imaginary part.
How to fix that in your code is more difficult, because you are the one who derived those equations. I would focus on those terms in your model where a number is raised to a non-integer power. First try to understand why the base of that exponent is negative. Why is the exponent itself some completely general fraction. You may need to deal with constraints on your model.

Iniciar sesión para comentar.


Sulaymon Eshkabilov
Sulaymon Eshkabilov el 2 de Sept. de 2024
Check your formulations for xh and xj
  1 comentario
Shaf
Shaf el 2 de Sept. de 2024
Hi, @Sulaymon Eshkabilov, thank you for the suggestion to recheck the formulation for xh and xj​. I found a small mistake and corrected it as follows:
xh = xa-xb+xc+xd+xe+xf+xg;
xj = ((8*(g(1)+hr+c*(hs+j))*(-(a*(-1+b)*g(1)*x(1)^2*y)/(L+g(1)*L*y))^(1/(1-b)))/(-1+b)+(2^(1/(-1+b))*l*pr*(2*L*(2*M+(-3+b))*(-(a*(-1+b)*g(1)*x(1)^2*y)/(L+g(1)*L*y))^(1/(1-b)) + 2*(-3+b)*L*(-(a*(-1+b)*g(1)*x(1)^2*y)/(L+g(1)*L*y))^(1/(1-b))))/((-1+b)*L*(1+g(1)*y))+x(1)*(2^(-2+b/(-1+b))*g(1)*y*((a*(1-b)*g(1)*x(1)^2*y)/(L+g(1)*L*y))^(b/(1-b)))/(L+g(1)*L*y)+x(1)^2*(2^(1/(-1+b))*l*pr*g(1)*y*(-3+b)*((a*(1-b)*g(1)*x(1)^2*y)/(L+g(1)*L*y))^(b/(1-b)))/((-1+b)*L*(1+g(1)*y)))/(2*x(1));
xk = (4*(c*Sr+Sv)/n^2+4*(S+c*T));
xl = (xk+x(1)*xh);
x(2) = sqrt(xl/-xj);
And for xh(i) and xj(i) :
xh(i)=(xa(i)-xb(i)+xc(i)+xd(i)+xe(i)+xf(i)+xg(i));
xj(i)=((8*(g(i-1)+hr+c*(hs+j))*(-(a*(-1+b)*g(i-1)*x(i-1)^2*y)/(L+g(i-1)*L*y))^(1/(1-b)))/(-1+b)+(2^(1/(-1+b))*l*pr*(2*L*(2*M+(-3+b))*(-(a*(-1+b)*g(i-1)*x(i-1)^2*y)/(L+g(i-1)*L*y))^(1/(1-b)) + 2*(-3+b)*L*(-(a*(-1+b)*g(i-1)*x(i-1)^2*y)/(L+g(i-1)*L*y))^(1/(1-b))))/((-1+b)*L*(1+g(i-1)*y))+x(i-1)*(2^(-2+b/(-1+b))*g(i-1)*y*((a*(1-b)*g(i-1)*x(i-1)^2*y)/(L+g(i-1)*L*y))^(b/(1-b)))/(L+g(i-1)*L*y)+x(i-1)^2*(2^(1/(-1+b))*l*pr*g(i-1)*y*(-3+b)*((a*(1-b)*g(i-1)*x(i-1)^2*y)/(L+g(i-1)*L*y))^(b/(1-b)))/((-1+b)*L*(1+g(i-1)*y)))/(2*x(i-1));
xk(i)=(4*(c*Sr+Sv)/n^2+4*(S+c*T));
xl(i)=(xk(i)+x(i-1)*xh(i));
x(i) =sqrt(xl(i)/-xj(i));
However, after running the corrected code, the result still contains an imaginary component.
Is there anything specific I should check in this formulation to identify the source of the imaginary part?

Iniciar sesión para comentar.

Categorías

Más información sobre Verification, Validation, and Test en Help Center y File Exchange.

Etiquetas

Productos


Versión

R2022a

Community Treasure Hunt

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

Start Hunting!

Translated by