\$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 / ball_rad 'ball radius' / 0.5 / ; midI = 1 + M / 2; midJ = (card(J) + 1) / 2; d(I,J) = sqrt(sqr((ord(I)-midI)/M) + sqr((ord(J)-midJ)/N)); ubnd(I,J) = [d(I,J) <= ball_rad]; v.up(ubnd(I,J)) = cheight - sqrt(ball_rad-sqr(d(I,J))); * 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), GMSPYTHONHOME gets set automatically * if pySetup=1 (default). If the %gams.sysdir%GMSPython directory does not exist on * these platforms pySetup will be reset to 0 and GAMS relies on an external Python * setup. \$if not setEnv GMSPYTHONHOME \$if %gams.pySetup%==1 \$abort.noError Embedded code Python not ready to be used \$if not setEnv GMSPYTHONHOME \$log --- Using external Python \$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: mi.sync_db[sym].copy_symbol(gams.db[sym]) 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; );