Imaginary Component in Total Profit Calculation

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!

Answer 1

[Shashi Kiran]

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.

[Shaf]

Hi @Shashi Kiran, thank you so much for your insights on my code. I'm not sure how to adjust the implementation to avoid this.

Should I be changing the parameter values to prevent the term from becoming negative, or is it more about rechecking and possibly revising the equation itself?

[Shashi Kiran]

@Shaf,

I recommend checking if there are any errors in your implementation of the equations.

If you could share the equations here, I can help you verify whether there are any errors in the implementation.

[Shaf]

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

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

Check your formulations for xh and xj

[Shaf]

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?