obstacle.gms : An Obstacle Problem

Description

```Program to calculate the position of a membrane pushed up through a
(rectangular) hole in a rigid plate; in addition, there are rigid
obstacle(s) inside the hole (perhaps not at the same level as the plate).
These obstacles can constrain the membrane from above or from below.

The correct position of the membrane (the position of minimum energy) is
determined by the minimization of a quadratic function of the membrane
position, subject to the constraints imposed by the hole and the obstacle.
The MCP below arises from the optimality conditions for this QP.

The MCP is solved for multiple force constants c using GUSS and
GAMSModelInstance. This is the first model in the library of this
model type that is solved this way.
```

Small Model of Type : MCP

Category : GAMS Model library

Main file : obstacle.gms

``````\$TITLE An Obstacle Problem

\$ontext
Program to calculate the position of a membrane pushed up through a
(rectangular) hole in a rigid plate; in addition, there are rigid
obstacle(s) inside the hole (perhaps not at the same level as the plate).
These obstacles can constrain the membrane from above or from below.

The correct position of the membrane (the position of minimum energy) is
determined by the minimization of a quadratic function of the membrane
position, subject to the constraints imposed by the hole and the obstacle.
The MCP below arises from the optimality conditions for this QP.

The MCP is solved for multiple force constants c using GUSS and
GAMSModelInstance. This is the first model in the library of this
model type that is solved this way.

Ciarlet, P G, The Finite Element Method for Elliptic Problems.
Elsevier Science, Studies in Mathematics and its Applications, 1978.
\$offtext

Sets
I       / 0 * 11 /
J       / 0 * 11 /
boundary(I,J)
interior(I,J)
ubnd(I,J) 'points with an upper bound';

boundary(I,J)\$[(ord(J) eq card(J)) or (ord(J) eq 1) or
(ord(I) eq card(I)) or (ord(I) eq 1)] = yes;
interior(I,J) = not boundary(I,J);

Scalars
xlo     / 0 /
xhi     / 1 /
ylo     / 0 /
yhi     / 1 /
dx, dy
c   'force constant' / 1 /
M   'number of divisions in I'
N   'number of divisions in J';

M = (card(I)-1);
N = (card(J)-1);
dx = (xhi - xlo) / M;
dy = (yhi - ylo) / N;

Parameters
d(I,J)   'distance of point from center of rectangle'
tmp(I,J)
ub(I,J);

Variable v(I,J)  'height of membrane';

Equation dv(I,J);

dv(I,J) ..
(dy/dx) * (2*v(I,J) - v(I+1,J) - v(I-1,J))
+ (dx/dy) * (2*v(I,J) - v(I,J+1) - v(I,J-1))
=N=
c * dx * dy;

Model obstacle / dv.v /;

* problem with a ball constraining from above
Scalar
midI, midJ
cheight  'height of center of ball' / .55 /
;
midI = 1 + M / 2;
midJ = (card(J) + 1) / 2;
d(I,J) = sqrt(sqr((ord(I)-midI)/M) + sqr((ord(J)-midJ)/N));

* the position v is fixed at zero on the boundary
v.fx(boundary) = 0;

v.l(interior(I,J)) = min(cheight, v.up(I,J));

Set s 'force scenarios' / s1*s5 /;
Parameter cval(s);

cval(s) = (ord(s)-1)*(1.5-0.5)/(card(s)-1) + 0.5;
display cval;

Set mattrib / system.GUSSModelAttributes /;

Parameter
vs(s,I,J)
srep(s, mattrib)  'model attibutes like modelstat etc'
o(*)              'GUSS options' / SkipBaseCase 1 /;

Set obstdict / s.scenario.'', o.opt.srep, c.param.cval, v.level.vs /;

solve obstacle using mcp scenario obstdict;

option vs:4:1:1; display vs;

* Now do the same with a GAMSModelInstance
*
* On the major platforms (Windows, Linux, Mac), GMSPYTHONLIB gets automatically set
* to use the internal Python installation in sysdir/GMSPython.
\$if not setEnv GMSPYTHONLIB \$abort.noError Embedded code Python not ready to be used
\$log --- Using Python library %sysEnv.GMSPYTHONLIB%

\$set solverlog
\$if set useSolverLog \$set solverlog output=sys.stdout

embeddedCode Python:
def solveMI(mi, symIn=[], symOut=[]):
for sym in symIn:
gams.db[sym].copy_symbol(mi.sync_db[sym])
mi.solve(%solverlog%)
for sym in symOut:
try:
gams.db[sym].clear() # Explicitly clear the symbol to ensure setting "writtenTo" flag for sym
mi.sync_db[sym].copy_symbol(gams.db[sym])
except:
pass
gams.wsWorkingDir = "."
pauseEmbeddedCode
abort\$execerror 'Python error. Check the log';

\$libInclude pyEmbMI miObst 'obstacle using mcp' -all_model_types=path c.Zero

option v:4:1:1;
loop(s,
c = cval(s);
continueEmbeddedCode:
gams.db['c'].copy_symbol(miObst.sync_db['c'])
solveMI(miObst,['c'],['v'])
pauseEmbeddedCode v
display v.l;
);
``````
GAMS Development Corp.
GAMS Software GmbH

General Information and Sales
U.S. (+1) 202 342-0180
Europe: (+49) 221 949-9170