$title Educational embedded complementarity system model (FERRIS43,SEQ=24)

$onText

  Embedded Complementarity System
  -------------------------------

  The problem to solve is:
    min_x  f(x,y)  st  g(x,y) <= 0
  plus the constraint that
    H(x,y,lambda) = 0
  where lambda is the multiplier on the g(x,y) <= 0 constraint.

  We use EMP to specify that the variable lambda appearing in the
  algebra is not a decision variable in any agent's model, but rather
  it is the dual variable to equation g.  N.B.: the equation g should
  be a typical constraint owned by an optimizing agent.

  We can describe the optimization model in two ways.  The "top-down"
  approach starts with a long-form solve statement and one minimizing
  agent owning all vars and equs.  The dualequ keyword removes H and y
  from this agent.  The "bottom-up" approach starts with an
  equilibrium system, adds a minimizing agent owning just the right
  equations and vars, and includes H and y as part of a separate VI
  agent.

  dualequ H y
  dualvar lambda g

References:
Ferris et al, An extended mathematical programming framework,
Computers and Chemical Engineering 33, p.1973-1982, 2009

Contributors: Jan-H. Jagla, November 2009.  Steve, 2017

$offText

variables obj, x, y;
positive variable lambda;
equations defobj, g, H;

defobj.. obj =e= sqr(x-y);
g.. y =g= x + 1;

H.. y + lambda =e= 2;

model ecs / defobj, g, H /;
file empinfo / '%emp.info%' /;

*-----------------------------------------------------------------------------
* Use long-form solve statement and peel off H and y via the dualequ directive
putclose empinfo
'dualequ H y' /
'dualvar lambda g' /
;
solve ecs using emp minimizing obj;

*------------------------------------------------------------------------
* Use short-form solve statement and build the equilibrium model in EMP,
* so H and y never get into minimizing agent's model
putclose empinfo
'equilibrium' /
' min obj x g defobj' /
' vi H y' /
' dualvar lambda g' /
;
solve ecs using EMP;

*------------------------------------------------------
* Write this model down manually and verify solution
equation dLdx;

dLdx.. ( - 2*(x - y))/(-1) + lambda =N= 0;

model mcp / g.lambda,dLdx.x,H.y /;
mcp.iterlim = 0;
solve mcp using MCP;
abort$(mcp.objval > 1e-6) 'Solutions not the same';