Imaginary Component in Total Profit Calculation

Question

Shaf 2024-9-2

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2149479-imaginary-component-in-total-profit-calculation

编辑： John D'Errico 2024-9-2，11:39

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!

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

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

Answer 1

Shashi Kiran 2024-9-2

1
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2149479-imaginary-component-in-total-profit-calculation#answer_1508679

在 MATLAB Online 中打开

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 个评论
显示 2更早的评论隐藏 2更早的评论

Shaf 2024-9-2，9:12

@Shashi Kiran

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 2024-9-2，11:38

编辑：John D'Errico 2024-9-2，11:39

在 MATLAB Online 中打开

@Shaf,

@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.

请先登录，再进行评论。

Answer 2

Sulaymon Eshkabilov 2024-9-2

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2149479-imaginary-component-in-total-profit-calculation#answer_1508644

Check your formulations for xh and xj

1 个评论
显示 -1更早的评论隐藏 -1更早的评论

Shaf 2024-9-2

在 MATLAB Online 中打开

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?

请先登录，再进行评论。