Question: μ₁ and μ₂ conditions in constrained optimization

I am optimizing an objective function subject to two constraints. I have obtained the solution, but I am unsure how to determine the conditions for the Lagrange multipliers μ1​ and μ2​. For example, using the KKT conditions, if μ1>0 then the solution is X, and if μ2>0 then the solution is Y. Could you guide me on how to find μ1 and μ2 along with their corresponding conditions?

restart

with(Optimization); with(plots); with(LinearAlgebra)

_local(Pi)

Pi

(1)
 

`π_m` := proc (i1) options operator, arrow; (w-i1)*(1/2+(1/2)*(i1-i2)/tau)*(1-(Pn-Pr)/(1-delta))-C0-(1/2)*Cr*(1/2+(1/2)*(i1-i2)/tau)^2*(1-(Pn-Pr)/(1-delta))^2+Ce*rho0*(1-(Pn-Pr)/(1-delta)) end proc

C1 := (2*rho0-1)*tau+i2 <= i1

(2*rho0-1)*tau+i2 <= i1

(2)

C2 := i1 <= tau+i2

i1 <= tau+i2

(3)

NULL

NULL

``

NULL

# No equality constraints
#
# Inequality constraints must be of the form f[i](i1) <= 0
#
# Thus
#
# (1) For C1
#     so
      f[1] := (i1) -> (2*rho0 - 1)*tau + i2-i1:
      
#
# (2) C2
#     so
      f[2] := (i1) -> i1-(tau + i2):
      
#

# Lagrangian (we want to maximize `&pi;_m` so to minimize -`&pi;_m`

L := -`&pi;_m`(i1) + add(f[i](i1)*mu[i], i=1..2):

dLdi1 := collect(diff(L, i1), [i1]):

KKT_conditions := [
                    seq(mu[i] >= 0, i=1..2),         # Dual feasibility conditions
                    dLdi1 = 0,                       # Stationarity condition
                    seq(``(f[i](i1)) <= 0, i=1..2),  # Primal feasibility conditions
                    add(mu[i]*f[i](i1) = 0, i=1..2)  # Complementary slackness
                  ]:

print~(KKT_conditions);

0 <= mu[1]

 

0 <= mu[2]

 

((1-(Pn-Pr)/(1-delta))/tau+(1/4)*Cr*(1-(Pn-Pr)/(1-delta))^2/tau^2)*i1+(1/2-(1/2)*i2/tau)*(1-(Pn-Pr)/(1-delta))-(1/2)*w*(1-(Pn-Pr)/(1-delta))/tau+(1/2)*Cr*(1/2-(1/2)*i2/tau)*(1-(Pn-Pr)/(1-delta))^2/tau-mu[1]+mu[2] = 0

 

``((2*rho0-1)*tau+i2-i1) <= 0

 

``(i1-tau-i2) <= 0

 

((2*rho0-1)*tau+i2-i1)*mu[1]+(i1-tau-i2)*mu[2] = 0

 

[]

(4)


Analysis of dLdp1

with(LargeExpressions):

DLDi1 := collect(dLdi1, i1, Veil[K]);

(1/4)*K[1]*i1-(1/4)*K[2]

(5)

# Thus dLdi1 = 0 verifies

isolate((5), i1)
 

i1 = K[2]/K[1]

(6)

# Where K1 and K2 are

i := 'i':
KS := [ seq(K[i] = Unveil[K](K[i]), i=1..LastUsed[K] ) ];

[K[1] = (-1+delta+Pn-Pr)*(Cr*Pn-Cr*Pr+Cr*delta+4*delta*tau-Cr-4*tau)/(tau^2*(-1+delta)^2), K[2] = (4*delta^2*tau^2*mu[1]-4*delta^2*tau^2*mu[2]+Cr*Pn^2*i2-Cr*Pn^2*tau-2*Cr*Pn*Pr*i2+2*Cr*Pn*Pr*tau+2*Cr*Pn*delta*i2-2*Cr*Pn*delta*tau+Cr*Pr^2*i2-Cr*Pr^2*tau-2*Cr*Pr*delta*i2+2*Cr*Pr*delta*tau+Cr*delta^2*i2-Cr*delta^2*tau+2*Pn*delta*i2*tau-2*Pn*delta*tau^2+2*Pn*delta*tau*w-2*Pr*delta*i2*tau+2*Pr*delta*tau^2-2*Pr*delta*tau*w+2*delta^2*i2*tau-2*delta^2*tau^2+2*delta^2*tau*w-8*delta*tau^2*mu[1]+8*delta*tau^2*mu[2]-2*Cr*Pn*i2+2*Cr*Pn*tau+2*Cr*Pr*i2-2*Cr*Pr*tau-2*Cr*delta*i2+2*Cr*delta*tau-2*Pn*i2*tau+2*Pn*tau^2-2*Pn*tau*w+2*Pr*i2*tau-2*Pr*tau^2+2*Pr*tau*w-4*delta*i2*tau+4*delta*tau^2-4*delta*tau*w+4*tau^2*mu[1]-4*tau^2*mu[2]+Cr*i2-Cr*tau+2*i2*tau-2*tau^2+2*tau*w)/(tau^2*(-1+delta)^2)]

(7)


Analysis of the two constraints

beta

beta

(8)

Cs := { seq(beta[i] = subs(i1=0, f[i](i1)), i=1..2) };

{beta[1] = (2*rho0-1)*tau+i2, beta[2] = -tau-i2}

(9)

 

# two constraints problem

Cmin_value := min(rhs~(Cs));  # the "true" constraints
Cmax_value := max(rhs~(Cs)):

Cmin_name  := min(lhs~(Cs));  # the "abstract" constraints
Cmax_name  := max(lhs~(Cs)):

Cmin_name  := LowerBound;  # and an even more abstract form
Cmax_name  := UpperBound:

Complementary_Slackness := mu[1]*(Cmin_name-i1) , mu[2]*(i1-Cmax_name);

L2     := -`&pi;_m`(i1) + add(Complementary_Slackness):
dL2di1 := diff(L2, i1):
dL2di1 := map(simplify, %, size);
 

min(-tau-i2, (2*rho0-1)*tau+i2)

 

min(beta[1], beta[2])

 

LowerBound

 

mu[1]*(LowerBound-i1), mu[2]*(i1-UpperBound)

 

(1/2)*(tau+i1-i2)*(-1+delta+Pn-Pr)/(tau*(-1+delta))+(1/2)*(-w+i1)*(-1+delta+Pn-Pr)/(tau*(-1+delta))+(1/4)*Cr*(tau+i1-i2)*(-1+delta+Pn-Pr)^2/(tau^2*(-1+delta)^2)-mu[1]+mu[2]

(10)

type(dL2di1, linear([i1, seq(mu[i], i=1..2)]));

SC_CS_sols := solve(
  {
    dL2di1 = 0,
    op([Complementary_Slackness] =~ 0)
  }
  ,
  {i1, seq(mu[i], i=1..2)}
):

true

 

Main: Entering solver with 3 equations in 3 variables
Main: attempting to solve as a linear system
Main: attempting to solve as a polynomial system
Main: Polynomial solver successful. Exiting solver returning 1 solution

 

# How many solutions did we get?

numelems([SC_CS_sols]);

# And those solutions are charecterized by

map(s -> if eval(lambda[1], s) = 0 then
           if eval(lambda[2], s) = 0 then
             "i1 belongs to interval (LowerBound, UpperBound)"
           else
             "i1 is equal to the UpperBound"
           end if
         else
           "i1 is equal to the LowerBound"
         end if
         , [SC_CS_sols]);

# So we get as expected the three possible situations.

3

 

["i1 is equal to the LowerBound", "i1 is equal to the LowerBound", "i1 is equal to the LowerBound"]

(11)

Cmin_value

min(-tau-i2, (2*rho0-1)*tau+i2)

(12)
 

 

Download I_1_Optimum_condition.mw

Please Wait...