$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';