# 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);

with(LargeExpressions):

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

# Thus dLdi1 = 0 verifies

isolate((5), i1)
 

# Where K₁ and K₂ are

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

beta

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

# 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);
 

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)}
):

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.

Cmin_value