Andiguys

65 Reputation

5 Badges

1 years, 79 days

MaplePrimes Activity


These are replies submitted by Andiguys

@acer 
Showing error "Error, (in assuming) when calling 'solve'. Received: 'cannot solve for an unknown function with other operations in its arguments'
"

Attaching sheet:

restart

kernelopts(version)

`Maple 2019.1, X86 64 WINDOWS, May 21 2019, Build ID 1399874`

(1)

B := w+t+((3*t+3*Pc-3*Pu)*upsilon+((-t-p1)*beta+a)*`ϕ`)/(3*eta); A := (1/3)*`ϕ`*(-beta*(p1+t)+a)-upsilon*Pu

K1 := (A(-2*U*upsilon+1)-2*U*t*upsilon^2+(-B+2*t)*upsilon)/(2*upsilon*(U*upsilon-1))

(1/2)*((1/3)*varphi(-2*U*upsilon+1)*(-beta(-2*U*upsilon+1)*(p1(-2*U*upsilon+1)+t(-2*U*upsilon+1))+a(-2*U*upsilon+1))-upsilon(-2*U*upsilon+1)*Pu(-2*U*upsilon+1)-2*U*t*upsilon^2+(-w+t-(1/3)*((3*t+3*Pc-3*Pu)*upsilon+((-t-p1)*beta+a)*varphi)/eta)*upsilon)/(upsilon*(U*upsilon-1))

(2)

K2 := simplify(solve(K1 = Pc, Pc))

(((-t(-2*U*upsilon+1)-p1(-2*U*upsilon+1))*beta(-2*U*upsilon+1)+a(-2*U*upsilon+1))*eta*varphi(-2*U*upsilon+1)-3*upsilon(-2*U*upsilon+1)*Pu(-2*U*upsilon+1)*eta-6*((U*t*eta+(1/2)*t-(1/2)*Pu)*upsilon+(-(1/2)*t+(1/2)*w)*eta+(1/6)*varphi*(-beta*p1-beta*t+a))*upsilon)/((6*U*eta+3)*upsilon^2-6*upsilon*eta)

(3)

``

ineq := -A/upsilon-t <= K2 and K2 <= (`&varphi;`*(-beta*p1+a)-A)/upsilon-t

 

You gave us the (new, updated) information that all variables
are positive. Literally, that means p1 as well (unless you now exclude it).

extra := indets(ineq,And(name,Not(constant))) >~ 0;

{0 < Pu, 0 < U, 0 < a, 0 < beta, 0 < eta, 0 < p1, 0 < t, 0 < upsilon, 0 < varphi, 0 < w}

(4)

 

ineq1 := op(1,ineq):

R1 := op(solve({ineq1},p1)) assuming extra[];

Error, (in assuming) when calling 'solve'. Received: 'cannot solve for an unknown function with other operations in its arguments'

 

ineq2 := op(2,ineq):

R2 := op(solve({ineq2},p1) assuming extra[]);

Error, (in assuming) when calling 'solve'. Received: 'cannot solve for an unknown function with other operations in its arguments'

 

Maybe this is all you are after...

`and`(R1, R2, p1>0);

R1 and R2 and 0 < p1

(5)

Or maybe something like this...

p1>0 and p1 <= min(rhs(R1),rhs(R2));

Error, invalid input: rhs received R1, which is not valid for its 1st argument, expr

 

Or maybe something like this...
simplify(solve(`and`(R1, R2, p1>0), p1)) assuming extra[];

Error, (in assuming) when calling 'simplify'. Received: 'invalid input: simplify uses a 1st argument, s, which is missing'

 
 

``

Download Q_isolate_acc_2_me.mw

@mmcdara Hi...if possible can you verify that the substitutions i made and the interpretation on optimal p1 which i have taken in above sheet is right or not? Thankyou

@mmcdara 

I have two questions:

  1. Is my approach for finding the optimal value of p1—based on whether the solution lies at the lower bound, upper bound, or within the interior—correct as per the syntax you provided?

  2. Are the lower and upper bounds that I substituted at the end of the sheet accurate? If not, what should they be?

Question_Lower_upper.mw

Status:

  • I've already formulated the KKT conditions.

  • Now I need help in solving for p1∗ in each case, depending on which constraints are binding.

Could someone help solve for p1∗​ under different combinations of active/passive constraints?

Looking for:

  • Analytical expressions for p1∗

  • Feasibility ranges per case

  • Guidance on structuring this systematically

@mmcdara 
Can you help with below code so that solve function execute fast?

restart

kernelopts(version)

`Maple 2019.1, X86 64 WINDOWS, May 21 2019, Build ID 1399874`

(1)

A := (1/3)*`&varphi;`*(-beta*p1+a)-upsilon*Pu

``

ineq := -((1/3)*`&varphi;`*(-beta*p1+a)-upsilon*Pu)/upsilon <= (2*U*beta*eta*nu*p1*`&varphi;`+6*Pu*U*eta*nu*upsilon-2*U*a*eta*nu*`&varphi;`+beta*eta*p1*`&varphi;`-beta*nu*p1*`&varphi;`+3*Pu*eta*upsilon-3*Pu*nu*upsilon-a*eta*`&varphi;`+a*nu*`&varphi;`+3*eta*nu*w)/(6*U*eta*nu^2+6*eta*nu-6*nu*upsilon+3*upsilon^2) and (2*U*beta*eta*nu*p1*`&varphi;`+6*Pu*U*eta*nu*upsilon-2*U*a*eta*nu*`&varphi;`+beta*eta*p1*`&varphi;`-beta*nu*p1*`&varphi;`+3*Pu*eta*upsilon-3*Pu*nu*upsilon-a*eta*`&varphi;`+a*nu*`&varphi;`+3*eta*nu*w)/(6*U*eta*nu^2+6*eta*nu-6*nu*upsilon+3*upsilon^2) <= (`&varphi;`*(-beta*p1+a)-A)/upsilon

 

extra := indets(ineq,And(name,Not(constant))) >~ 0;

{0 < Pu, 0 < U, 0 < a, 0 < beta, 0 < eta, 0 < nu, 0 < p1, 0 < upsilon, 0 < varphi, 0 < w}

(2)

 

ineq1 := op(1,ineq):

R1 := op(solve({ineq1},p1)) assuming extra[];

ineq2 := op(2,ineq):

R2 := op(solve({ineq2},p1) assuming extra[]);



`and`(R1, R2, p1>0);



p1>0 and p1 <= min(rhs(R1),rhs(R2));


simplify(solve(`and`(R1, R2, p1>0), p1)) assuming extra[];

 

``

Download Q_isolate_acc_2.mw

@acer 

Apologies for the confusion earlier — I should have mentioned that all variables are positive and should have included both constraints together from the beginning.

Initially, my intent was to understand the correct syntax, and once I had that working, I tried applying it to a slightly different objective function. However, with this new function, the same syntax is now taking a very long time to solve, and I haven’t received any result yet.

Would you be able to suggest what changes or simplifications I can make to improve the computation time for this case?

@acer Taking too much time to solve? Why?

restart

kernelopts(version)

`Maple 2019.1, X86 64 WINDOWS, May 21 2019, Build ID 1399874`

(1)

B := w+Ce; A := (1/3)*`&varphi;`*(-beta*p1+a)-upsilon*Pu

ineq := (2*U*beta*eta*nu*p1*`&varphi;`+6*Pu*U*eta*nu*upsilon-2*U*a*eta*nu*`&varphi;`+beta*eta*p1*`&varphi;`-beta*nu*p1*`&varphi;`+3*Pu*eta*upsilon-3*Pu*nu*upsilon-a*eta*`&varphi;`+a*nu*`&varphi;`+3*eta*nu*w)/(6*U*eta*nu^2+6*eta*nu-6*nu*upsilon+3*upsilon^2) <= (2*`&varphi;`*(-beta*p1+a)*(1/3)+upsilon*Pu)/upsilon

(2*U*beta*eta*nu*p1*varphi+6*Pu*U*eta*nu*upsilon-2*U*a*eta*nu*varphi+beta*eta*p1*varphi-beta*nu*p1*varphi+3*Pu*eta*upsilon-3*Pu*nu*upsilon-a*eta*varphi+a*nu*varphi+3*eta*nu*w)/(6*U*eta*nu^2+6*eta*nu-6*nu*upsilon+3*upsilon^2) <= ((2/3)*varphi*(-beta*p1+a)+upsilon*Pu)/upsilon

(2)

temp := collect(ineq,p1);

(2*U*beta*eta*nu*varphi+beta*eta*varphi-beta*nu*varphi)*p1/(6*U*eta*nu^2+6*eta*nu-6*nu*upsilon+3*upsilon^2)+(6*Pu*U*eta*nu*upsilon-2*U*a*eta*nu*varphi+3*Pu*eta*upsilon-3*Pu*nu*upsilon-a*eta*varphi+a*nu*varphi+3*eta*nu*w)/(6*U*eta*nu^2+6*eta*nu-6*nu*upsilon+3*upsilon^2) <= -(2/3)*varphi*beta*p1/upsilon+((2/3)*varphi*a+upsilon*Pu)/upsilon

(3)

# lhs and rhs are both of type `+`, so select/remove make sense
new := simplify(temp - remove(depends,lhs(temp),p1) - select(depends,rhs(temp),p1));

(2/3)*p1*varphi*((1/2)*upsilon^2+((1/2)*U*nu*eta-(5/4)*nu+(1/4)*eta)*upsilon+nu*eta*(U*nu+1))*beta/(upsilon*((1/2)*upsilon^2-nu*upsilon+nu*eta*(U*nu+1))) <= (1/6)*(3*Pu*upsilon^3+((-6*Pu*U*eta-3*Pu)*nu+2*varphi*a-3*Pu*eta)*upsilon^2+(6*U*Pu*nu^2*eta+((2*U*a*varphi+6*Pu-3*w)*eta-5*varphi*a)*nu+a*eta*varphi)*upsilon+4*a*nu*eta*varphi*(U*nu+1))/(upsilon*((1/2)*upsilon^2-nu*upsilon+nu*eta*(U*nu+1)))

(4)

#Taking too much time to solve???
simplify(solve(new,p1));

Download Q_isolate_acc.mw

@acer Thankyou

@acer Thank you.

If I replace the condition a ≤ b with a range like c ≤ a ≤ b, where a, b, and c are functions of p1, how can I isolate p1 to determine its range?

Note: provided all parameters are positive.

restart

kernelopts(version)

`Maple 2019.1, X86 64 WINDOWS, May 21 2019, Build ID 1399874`

(1)

B := w+Ce; A := (1/3)*`&varphi;`*(-beta*p1+a)-upsilon*Pu

ineq := -A/upsilon <= (B-A*(2*U*upsilon+1)/upsilon)/(2*(U*upsilon+1)) and (B-A*(2*U*upsilon+1)/upsilon)/(2*(U*upsilon+1)) <= (`&varphi;`*(-beta*p1+a)-A)/upsilon

-((1/3)*varphi*(-beta*p1+a)-upsilon*Pu)/upsilon <= (w+Ce-((1/3)*varphi*(-beta*p1+a)-upsilon*Pu)*(2*U*upsilon+1)/upsilon)/(2*U*upsilon+2) and (w+Ce-((1/3)*varphi*(-beta*p1+a)-upsilon*Pu)*(2*U*upsilon+1)/upsilon)/(2*U*upsilon+2) <= ((2/3)*varphi*(-beta*p1+a)+upsilon*Pu)/upsilon

(2)

temp := collect(ineq,p1);

(1/3)*varphi*beta*p1/upsilon-((1/3)*varphi*a-upsilon*Pu)/upsilon <= (1/3)*varphi*beta*(2*U*upsilon+1)*p1/(upsilon*(2*U*upsilon+2))+(w+Ce-((1/3)*varphi*a-upsilon*Pu)*(2*U*upsilon+1)/upsilon)/(2*U*upsilon+2) and (1/3)*varphi*beta*(2*U*upsilon+1)*p1/(upsilon*(2*U*upsilon+2))+(w+Ce-((1/3)*varphi*a-upsilon*Pu)*(2*U*upsilon+1)/upsilon)/(2*U*upsilon+2) <= -(2/3)*varphi*beta*p1/upsilon+((2/3)*varphi*a+upsilon*Pu)/upsilon

(3)

# lhs and rhs are both of type `+`, so select/remove make sense
new := simplify(temp - remove(depends,lhs(temp),p1) - select(depends,rhs(temp),p1));

Error, invalid input: lhs received (1/3)*varphi*beta*p1/upsilon-((1/3)*varphi*a-upsilon*Pu)/upsilon <= (1/3)*varphi*beta*(2*U*upsilon+1)*p1/(upsilon*(2*U*upsilon+2))+(w+Ce-((1/3)*varphi*a-upsilon*Pu)*(2*U*upsilon+1)/upsilon)/(2*U*upsilon+2) and (1/3)*varphi*beta*(2*U*upsilon+1)*p1/(upsilon*(2*U*upsilon+2))+(w+Ce-((1/3)*varphi*a-upsilon*Pu)*(2*U*upsilon+1)/upsilon)/(2*U*upsilon+2) <= -(2/3)*varphi*beta*p1/upsilon+((2/3)*varphi*a+upsilon*Pu)/upsilon, which is not valid for its 1st argument, expr

 

simplify(solve(new,p1));

Error, invalid input: simplify uses a 1st argument, s, which is missing

 
 

``

Download Q_isolate_acc.mw

@mmcdara So, from the sheet, I gathered that when 1×10^-10 ≤ U_0   and upsilon ≤ 1, both constraints are satisfied. Is that correct?

@vv In constraint p1 is in RHS and LHS, how to bring it to one side? Thats why optimization code not reading it.

    C1 := -0.3500000000 p1 <= -15289.70833 + 0.1837500000 p1

    C2 := 0.1837500000 p1 <= -0.4083333334 p1 + 16956.66667

Stil showing error : "Error, (in Optimization:-NLPSolve) no feasible point found for the linear constraints"

@mmcdara Similarly, I attempted to merge the two plots, but in the second case, the plot appears identical to the first one. Could you help me identify what might have gone wrong? I’ve attached the actual plot for Case 2, created on a separate sheet, where you can see that the value of Pc starts around 2820. However, when I combine both cases into a single sheet, the Pc values for Case 2 seem to be taken from Plot 1 instead.

Problem mentioned in sheet: Q_N1_mmcdara.mw

Case 2: 

@mmcdara I have two plots representing two different cases, both with TR on the x-axis and TP on the y-axis. I’d like to combine them into a single plot, using different colored lines to distinguish each case. What code should I use to do this?

restart

_local(D)

DATA := [Cm = 8000, U[0] = 5, Cr = 2000, a = 19000, U = 500, w = 2000, lambda = .9, tau0 = .4, k = 10000, Pu = 500, beta = .8, upsilon = .8, `&varphi;` = .8, t = 1000]

NULL

L_wout := proc (p1, p2) options operator, arrow; (p1-Cm)*(-beta*p1+a)+lambda*(((1/3)*`&varphi;`*(-beta*p1+a)-upsilon*(Pu-Pc))*(-epsilon*p1-Cr+p2-w)+(-beta*p2+a-(1/3)*`&varphi;`*(-beta*p1+a)+upsilon*(Pu-Pc))*(p2-Cm)+k*((1/3)*`&varphi;`*(-beta*p1+a)-upsilon*(Pu-Pc)-tau0*(-beta*p1+a))) end proc

R_wout := proc (Pc) options operator, arrow; lambda*(w*((1/3)*`&varphi;`*(-beta*p1+a)-upsilon*(Pu-Pc))+epsilon*p1*((1/3)*`&varphi;`*(-beta*p1+a)-upsilon*(Pu-Pc))-Pc*((1/3)*`&varphi;`*(-beta*p1+a)-upsilon*(Pu-Pc))-U*((1/3)*`&varphi;`*(-beta*p1+a)-upsilon*(Pu-Pc))+U[0]*((1/3)*`&varphi;`*(-beta*p1+a)-upsilon*(Pu-Pc))^2) end proc

C1 := {-((1/3)*`&varphi;`*(-beta*p1+a)-upsilon*Pu)/upsilon <= Pc}

C2 := {Pc <= (2*`&varphi;`*(-beta*p1+a)*(1/3)+upsilon*Pu)/upsilon}

C_11 := subs(p1 = (-6*beta*k*lambda*tau0*upsilon*U[0]+Cm*beta*lambda*`&varphi;`-6*Cm*beta*upsilon*U[0]-3*Cm*epsilon*lambda*upsilon-Cr*beta*lambda*`&varphi;`+3*Cr*epsilon*lambda*upsilon-3*Pu*epsilon*lambda*upsilon-3*U*epsilon*lambda*upsilon+a*epsilon*lambda*`&varphi;`-6*beta*k*lambda*tau0+beta*k*lambda*`&varphi;`-beta*lambda*w*`&varphi;`-3*epsilon*k*lambda*upsilon+6*epsilon*lambda*upsilon*w-6*a*upsilon*U[0]-6*Cm*beta-6*a)/(2*(beta*epsilon*lambda*`&varphi;`-3*epsilon^2*lambda*upsilon-6*beta*upsilon*U[0]-6*beta)), C1)

C_21 := subs(p1 = (-6*beta*k*lambda*tau0*upsilon*U[0]+Cm*beta*lambda*`&varphi;`-6*Cm*beta*upsilon*U[0]-3*Cm*epsilon*lambda*upsilon-Cr*beta*lambda*`&varphi;`+3*Cr*epsilon*lambda*upsilon-3*Pu*epsilon*lambda*upsilon-3*U*epsilon*lambda*upsilon+a*epsilon*lambda*`&varphi;`-6*beta*k*lambda*tau0+beta*k*lambda*`&varphi;`-beta*lambda*w*`&varphi;`-3*epsilon*k*lambda*upsilon+6*epsilon*lambda*upsilon*w-6*a*upsilon*U[0]-6*Cm*beta-6*a)/(2*(beta*epsilon*lambda*`&varphi;`-3*epsilon^2*lambda*upsilon-6*beta*upsilon*U[0]-6*beta)), C2)

C_1 := subs(DATA, C_11)

C_2 := subs(DATA, C_21)

NULL

TRC_1 := proc (Pc) options operator, arrow; subs(p1 = (-6*beta*k*lambda*tau0*upsilon*U[0]+Cm*beta*lambda*`&varphi;`-6*Cm*beta*upsilon*U[0]-3*Cm*epsilon*lambda*upsilon-Cr*beta*lambda*`&varphi;`+3*Cr*epsilon*lambda*upsilon-3*Pu*epsilon*lambda*upsilon-3*U*epsilon*lambda*upsilon+a*epsilon*lambda*`&varphi;`-6*beta*k*lambda*tau0+beta*k*lambda*`&varphi;`-beta*lambda*w*`&varphi;`-3*epsilon*k*lambda*upsilon+6*epsilon*lambda*upsilon*w-6*a*upsilon*U[0]-6*Cm*beta-6*a)/(2*beta*epsilon*lambda*`&varphi;`-6*epsilon^2*lambda*upsilon-12*beta*upsilon*U[0]-12*beta), R_wout(Pc)) end proc

TRC := subs(DATA, TRC_1(Pc))

``

NULL

L_w := proc (p1, p2) options operator, arrow; (p1+t-Cm)*(a-beta*(p1+t))+lambda*(((1/3)*`&varphi;`*(a-beta*(p1+t))-upsilon*(Pu-Pc-t))*(-epsilon*p1-Cr+p2-w)+(a-beta*(p2+t)-(1/3)*`&varphi;`*(a-beta*(p1+t))+upsilon*(Pu-Pc-t))*(p2+t-Cm)+k*((1/3)*`&varphi;`*(a-beta*(p1+t))-upsilon*(Pu-Pc-t)-tau0*(a-beta*(p1+t)))) end proc

R_w := proc (Pc) options operator, arrow; lambda*((w+t)*((1/3)*`&varphi;`*(a-beta*(p1+t))-upsilon*(Pu-Pc-t))+epsilon*p1*((1/3)*`&varphi;`*(a-beta*(p1+t))-upsilon*(Pu-Pc-t))-(Pc+t)*((1/3)*`&varphi;`*(a-beta*(p1+t))-upsilon*(Pu-Pc-t))-U*((1/3)*`&varphi;`*(a-beta*(p1+t))-upsilon*(Pu-Pc-t))+U[0]*((1/3)*`&varphi;`*(a-beta*(p1+t))-upsilon*(Pu-Pc-t))^2) end proc

NULL

D1 := {-``*((1/3)*`&varphi;`*(-beta*(p1+t)+a)-upsilon*Pu)/upsilon-t <= Pc}

D2 := {Pc <= (2*`&varphi;`*(-beta*(p1+t)+a)*(1/3)+upsilon*Pu)/upsilon-t}

D_11 := subs(p1 = (-6*beta*k*lambda*tau0*upsilon*U[0]-beta*epsilon*lambda*t*`&varphi;`+Cm*beta*lambda*`&varphi;`-6*Cm*beta*upsilon*U[0]-3*Cm*epsilon*lambda*upsilon-Cr*beta*lambda*`&varphi;`+3*Cr*epsilon*lambda*upsilon-3*Pu*epsilon*lambda*upsilon-3*U*epsilon*lambda*upsilon+a*epsilon*lambda*`&varphi;`-6*beta*k*lambda*tau0+beta*k*lambda*`&varphi;`-beta*lambda*t*`&varphi;`-beta*lambda*w*`&varphi;`+12*beta*t*upsilon*U[0]-3*epsilon*k*lambda*upsilon+6*epsilon*lambda*t*upsilon+6*epsilon*lambda*upsilon*w-6*a*upsilon*U[0]-6*Cm*beta+12*beta*t-6*a)/(2*(beta*epsilon*lambda*`&varphi;`-3*epsilon^2*lambda*upsilon-6*beta*upsilon*U[0]-6*beta)), D1)

D_21 := subs(p1 = (-6*beta*k*lambda*tau0*upsilon*U[0]-beta*epsilon*lambda*t*`&varphi;`+Cm*beta*lambda*`&varphi;`-6*Cm*beta*upsilon*U[0]-3*Cm*epsilon*lambda*upsilon-Cr*beta*lambda*`&varphi;`+3*Cr*epsilon*lambda*upsilon-3*Pu*epsilon*lambda*upsilon-3*U*epsilon*lambda*upsilon+a*epsilon*lambda*`&varphi;`-6*beta*k*lambda*tau0+beta*k*lambda*`&varphi;`-beta*lambda*t*`&varphi;`-beta*lambda*w*`&varphi;`+12*beta*t*upsilon*U[0]-3*epsilon*k*lambda*upsilon+6*epsilon*lambda*t*upsilon+6*epsilon*lambda*upsilon*w-6*a*upsilon*U[0]-6*Cm*beta+12*beta*t-6*a)/(2*(beta*epsilon*lambda*`&varphi;`-3*epsilon^2*lambda*upsilon-6*beta*upsilon*U[0]-6*beta)), D2)

NULL

NULL

D_1 := solve(subs(DATA, D_11), Pc)

D_2 := subs(DATA, D_21)

NULL

NULL

TRC_2 := proc (Pc) options operator, arrow; subs(p1 = (-6*beta*k*lambda*tau0*upsilon*U[0]-beta*epsilon*lambda*t*`&varphi;`+Cm*beta*lambda*`&varphi;`-6*Cm*beta*upsilon*U[0]-3*Cm*epsilon*lambda*upsilon-Cr*beta*lambda*`&varphi;`+3*Cr*epsilon*lambda*upsilon-3*Pu*epsilon*lambda*upsilon-3*U*epsilon*lambda*upsilon+a*epsilon*lambda*`&varphi;`-6*beta*k*lambda*tau0+beta*k*lambda*`&varphi;`-beta*lambda*t*`&varphi;`-beta*lambda*w*`&varphi;`+12*beta*t*upsilon*U[0]-3*epsilon*k*lambda*upsilon+6*epsilon*lambda*t*upsilon+6*epsilon*lambda*upsilon*w-6*a*upsilon*U[0]-6*Cm*beta+12*beta*t-6*a)/(2*beta*epsilon*lambda*`&varphi;`-6*epsilon^2*lambda*upsilon-12*beta*upsilon*U[0]-12*beta), R_w(Pc)) end proc

T_R := subs(DATA, TRC_2(Pc))

obj   := TRC:
cstr  := [C_1[], C_2[]]:
H     := lambda__1*(lhs(cstr[1])-Pc), lambda__2*(Pc-rhs(cstr[2])):
lag   := obj + add(H):
sols  := [solve([diff(lag, Pc), H], {Pc, lambda__1, lambda__2})]:
sols  := map(s -> if eval(lambda__1^2+lambda__2^2, s) = 0 then NULL else s end if, sols):
objs  := map(s -> eval(diff(obj, Pc), s), sols):
pairs := map(s -> [eval(Pc, s), eval(TRC, Pc=eval(Pc, s))], sols):

print():print():

TP_1   := eval(L_wout(p1, p2), sols[1]):
TP     := collect(subs(DATA, TP_1), [p1, p2]):
sol    := solve(diff~(TP, [p1, p2]), {p1, p2}):
triple := eval([p1, p2, TP], sol):

 

(1)

ob   := T_R:
cst  := [D_1[], D_2[]]:
H_1     := lambda__1*(lhs(cst[1])-Pc), lambda__2*(Pc-rhs(cst[2])):
lag_1   := ob + add(H_1):
sols_1  := [solve([diff(lag_1, Pc), H], {Pc, lambda__1, lambda__2})]:
sols_1  := map(s_1 -> if eval(lambda__1^2+lambda__2^2, s_1) = 0 then NULL else s_1 end if, sols_1):
objs_1  := map(s_1 -> eval(diff(ob, Pc), s_1), sols_1):
pairs_1 := map(s_1 -> [eval(Pc, s_1), eval(T_R, Pc=eval(Pc, s_1))], sols_1):

print():print():

TP_2   := eval(L_w(p1, p2), sols_1[1]):
TP_21     := collect(subs(DATA, TP_2), [p1, p2]):
sol_1    := solve(diff~(TP_21, [p1, p2]), {p1, p2}):
triple_1 := eval([p1, p2, TP_21], sol_1):

 

(2)

with(plots)

display(plot([pairs[1][2], triple[3], varepsilon = 0 .. .2], color = blue, linestyle = 3, labels = [TR, 'TP']), seq(textplot([(eval([pairs[1][2], triple[3]], varepsilon = val))[], nprintf("%1.04f", val)], align = {right}), val = 0 .. .2, 0.4e-1), seq(pointplot([eval([pairs[1][2], triple[3]], varepsilon = val)], symbol = solidcircle, symbolsize = 10, color = blue), val = 0 .. .2, 0.4e-1), size = [500, 400])

 

display(plot([pairs_1[1][2], triple_1[3], varepsilon = 0 .. .2], color = blue, linestyle = 3, labels = [TR, 'TP']), seq(textplot([(eval([pairs_1[1][2], triple_1[3]], varepsilon = val))[], nprintf("%1.04f", val)], align = {right}), val = 0 .. .2, 0.4e-1), seq(pointplot([eval([pairs_1[1][2], triple_1[3]], varepsilon = val)], symbol = solidcircle, symbolsize = 10, color = blue), val = 0 .. .2, 0.4e-1), size = [500, 400])

 
 

``

Download Q_N1.mw

I have a doubt regarding the result. When I run the optimization, I get a different value for TR compared to what's shown in the graph. For example, at varepsilon = 0.09, the graph indicates that TR is around 102,000, but when I solve it using NLP optimization, the TR value comes out to be approximately 6.03 × 10^7. Interestingly, the value of Pc is the same in both the graph and the optimized result. Could you help me understand what might have gone wrong? I’m also attaching a screenshot of the code I used, which I placed below the Case 1 scenario in your example.

@mmcdara 

Thank you for sharing the syntax. I have two questions:

  1. Just a Doubt: Are we optimizing the TR and TP functions? I couldn't find any maximize or NLPsolve  functions in the code you provided.

  2. How can I combine the graphs for Case 1 and Case 2 into a single plot? Specifically, I’d like to merge the plots for P1​, Pc​, and ε from both cases into one graph, and do the same for the TR and TP graphs. I’ve attached the spreadsheet for both cases.
    Case_1.mw
    Case_2.mw

1 2 3 4 5 6 7 Page 3 of 9