The purpose of this post is to detail work done by Dan Melanz and I on the topic of solving frictional contact problems using the General Algebraic Modeling System (GAMS)

The best way to describe GAMS is that it is a language and framework for writing mathematical programming and optimization problems. GAMS can solve many different classes of problems, in this post the Nonlinear Programming (NLP), Quadratically Constrained Program (QCP) and Extended Mathematical Program (EMP) problem types will be discussed.

## References

These papers provide more details about the various formulations used in this post.

#### Primal Model:

V. Acary and F. Cadoux, “Applications of an Existence Result for the Coulomb Friction Problem,” Recent Advances in Contact Mechanics, pp. 45–66, 2013.

V. Acary, F. Cadoux, C. Lemaréchal, and J. Malick, “A formulation of the linear discrete Coulomb friction problem via convex optimization,” Z. angew. Math. Mech., vol. 91, no. 2, pp. 155–175, Feb. 2011.

F. Bertails-Descoubes, F. Cadoux, G. Daviet, and V. Acary, “A nonsmooth Newton solver for capturing exact Coulomb friction in fiber assemblies,” ACM Transactions on Graphics, vol. 30, no. 1, pp. 1–14, Jan. 2011.

#### Dual Model:

M. Anitescu and G. D. Hart, “A constraint-stabilized time-stepping approach for rigid multibody dynamics with joints, contact and friction,” Int. J. Numer. Meth. Engng., vol. 60, no. 14, pp. 2335–2371, Aug. 2004.

A. Tasora and M. Anitescu, “A complementarity-based rolling friction model for rigid contacts,” Meccanica, vol. 48, no. 7, pp. 1643–1659, Sep. 2013.

## Equations of Motion (EOM)

#### Variables and symbols

• $M$ - Mass Matrix
• $D$ - Jacobian Matrix
• $q$ - Positions
• $v$ - Velocities
• $\widehat{\gamma}$ - Reaction Forces
• $\gamma$ - Reaction Impulses
• $h$ - Time Step
• $\mu$ - Friction Constant
• $N$ - Shur Complement Matrix
• $r$ - Shur Complement Right Hand Side
• $\Upsilon_i$ - Friction Code
• $\Upsilon_i^{\circ}$ - Polar Friction Cone

Below are the Equations for motion for a system of rigid bodies that have contact (unilateral) and joint (bilateral) constraints.
The second equation is a simple force balance stating $M \bf{\dot v} =f$ where the forces comprise of applied forces such as gravity, reaction forces from joints and frictional contact forces. The next equation states that the joints must be satisfied at the position level (No drift). The complementarity condition states that the lagrange multiplier $\widehat{\gamma}_{i,n}$ is positive when the gap $\Phi_{i}$ is zero and vice versa, a zero reaction force means that the gap is positive.

#### Force-Acceleration Form


#### Discretized Impulse-Velocity Form

Here the equations are similar, but the lagrange multipliers become reaction impulses rather than forces. Also stabilization terms are added for both the contacts and the joints that correct for penetrations. The relaxation term can be thought of as a “tilting” function that makes the lagrange multipliers perpendicular to the friction cone.


## Dynamics in three ways

I am going to discuss three different models that can be used for solving these equations within GAMS. The first is solving the equations of motion directly using an EMP formulation. The second will solve a quadratic optimization problem with the unknowns being our lagrange multipliers. The third will solve a similar optimization problem where the unknowns are the new velocities.

#### Setting up the basic GAMS model

I am going to assume that GAMS has been installed with a proper license and all of your solvers are in working order.

This is the basic framework that I will be using for all models. The main idea is to limit as much as possible the information written to disk, using SolveLink=5 keeps gams in memory and makes it run faster. This makes it useful when running GAMS in the loop with a dynamics engine. Here is common code that is needed for all models

Here the set constraints represents the total constraints in the system of equation. For example, one contact has three constraints, one normal and two tangential. The set contacts has a one to one mapping with the total number of contacts. The set dofs is a list containing all of the degrees of freedom for the system. Each body has 6 entries and they are organized in a linear fashion, every sixth entry represents a different body.

#### Solving The Quadratic Optimization Problem (Dual Form)

Here we transform the EOM to the following Cone Complementarity Problem (CCP)


This CCP leads to the following quadratic optimization problem with conic constraints


This model can be specified in GAMS in the following manner:

Using command line arguments the type of model can be changed between QCP and NLP. Different types of solvers can be used with each type of model.

#### Solving The Quadratic Optimization Problem (Primal Form)

This version of the quadratic optimization problem solves for the new velocities


This model is slightly different than the other two in that it uses an iterative scheme to update the relaxation term. We found that doing 10 iterations was usually enough to converge to a solution.

## Solvers

This table represents the ist of solvers that will be used for each model, the goal is to understand how they perform and if they provide different solutions

Solver Dual qcp Dual nlp Primal qcp Primal nlp EOM
Conopt YES YES YES YES NO
Cplexd YES NO YES NO NO
Gurobi YES NO YES NO NO
Ipopt YES YES YES YES NO
minos YES YES YES YES NO
knitro YES YES YES YES NO
lindo YES NO YES NO NO
pathnlp NO YES NO YES NO
miles NO YES NO NO YES
path NO YES NO NO YES

## Models

Here several models will be described and detals provided on the solutions provided by different solvers

#### Box on an inclined plane

Results for normal force and tangential friction force for an angle of 0 degrees. Mass of the block is 1.0 kg

Solver dual qcp dual nlp primal qcp primal nlp EOM emp
Conopt 9.81 9.81 9.81 9.81 0
Cplexd 9.699811889 0 9.81000067 0 0
Gurobi 10.00146358 0 9.812432398 0 0
Ipopt 9.81004453 9.81004453 9.810017486 9.810017486 0
minos 9.810051932 9.810051932 9.81 9.81 0
knitro 9.812031608 9.812031608 10.33712774 10.33712774 0
lindo 39.24048097 0 9.810070548 0 0
pathnlp 0 9.809653585 0 9.80965966 0
miles 0 0 0 0 9.81
path 0 0 0 0 9.809913388
Solver dual qcp dual nlp primal qcp primal nlp EOM emp
Conopt 0 0 1.63E-16 1.63E-16 0
Cplexd 9.15E-12 0 3.68E-20 0 0
Gurobi 4.82E-11 0 -2.35E-16 0 0
Ipopt -2.03E-16 -2.03E-16 -6.94E-07 -6.94E-07 0
minos 0 0 0 0 0
knitro -2.44E-16 -2.44E-16 -1.18E-17 -1.18E-17 0
lindo 6.36E-14 0 1.21E-19 0 0
pathnlp 0 0 0 0 0
miles 0 0 0 0 0
path 0 0 0 0 0

Results for normal force and tangential friction force for an angle of 25 degrees. Mass of the block is 1.0 kg

Solver dual qcp dual nlp primal qcp primal nlp EOM emp
Conopt 8.890880015 8.890880015 8.903193995 8.903193995 0
Cplexd 8.887833911 0 8.890880643 0 0
Gurobi 8.516675409 0 8.943110478 0 0
Ipopt 8.890897594 8.890897594 8.900874991 8.900874991 0
minos 8.890898214 8.890898214 8.891385779 8.891385779 0
knitro 8.890953799 8.890953799 9.988990071 9.988990071 0
lindo 35.56570555 0 8.891402252 0 0
pathnlp 0 8.890879976 0 8.890982302 0
miles 0 0 0 0 8.890880015
path 0 0 0 0 8.890880079
Solver dual qcp dual nlp primal qcp primal nlp EOM emp
Conopt 4.145881897 4.145881897 4.141172311 4.141172311 0
Cplexd 4.146100806 0 4.145881789 0 0
Gurobi 4.346501223 0 4.120027138 0 0
Ipopt 4.145873977 4.145873977 4.135891029 4.135891029 0
minos 4.146010481 4.146010481 4.145469386 4.145469386 0
knitro 4.145843537 4.145843537 3.742626642 3.742626642 0
lindo 16.5822553 0 4.145673861 0 0
pathnlp 0 4.145881917 0 4.145831445 0
miles 0 0 0 0 4.145881897
path 0 0 0 0 4.145881833

#### Timing Comparison

In this test there is a stack of spheres in a 4x20x4 configuration in contact with a flat plate. There are a total of 800 contacts, all times are in seconds. The Lindo solver did not finish after running for more than 300 seconds. s

Solver dual qcp dual nlp primal qcp primal nlp EOM emp
Conopt 0.232 0.231 136.702 137.543 0
Cplexd 0 0 0.1 0 0
Gurobi 0.24 0 0.16 0 0
Ipopt 1.127 1.116 0.385 0.409 0
minos 1.726 1.718 0.031 0.034 0
knitro 72.836 73.141 0.024 0.026 0
lindo DNF 0 0.077 0 0
pathnlp 0 0.049 0 0.021 0
miles 0 0 0 0 0.073
path 0 0 0 0 0.039