Nonlinear least squares and polygon area

I found this paper that describes preserving the global area of a polygon given some deformation (section 5): http://www.kunzhou.net/publications/2DShape.pdf

I'm trying to do something very similar. Right now I have a linear least squares solver with inequality constraints. eg: given a matrix $A$, find the Lagrangian/impulse ($\lambda$) to "solve" the system subject to some hard limits ($G\vec{\lambda} \le h$)

$min || A \vec{\lambda} - b ||^2$ s.t. $G\vec{\lambda} \le h$

This is a classic linear least squares with inequality constraints problem, as explained in this paper: ftp://reports.stanford.edu/pub/cstr/reports/cs/tr/69/134/CS-TR-69-134.pdf, and the implementation I have seems to work right now.

In my case $\vec{\lambda}$ is actually changes to the x,y components of the vertices of a polygon. eg:

$\vec{x} = [x_1, y_1, x_2, y_2, x_3, y_3, ... ]^T$

$\vec{x_{i+1}} = \vec{x_i} + \vec{\lambda}$

To what I have right now I'm trying to add additional hard limits that the area of the polygon has to be conserved. The paper I linked at the start has an algorithm for setting up and solving the area preserving constraint.

So I have a few questions:

1. I'm not sure how to combine the area preserving constraint (which is a nonlinear hard equality) with my existing linear inequality constraints in any meaningful way.
2. The paper uses an iterative Gauss-Newton solver, which if I understand things correctly is an interior point method. What I have right now uses Dantzig's basis exchange algorithm. Is there a way to adapt the basis exchange algorithm to work with the nonlinear constraints?
3. If I have to use an iterative solver as in that paper, do I realistically have to be concerned about divergence or anything weird like that? This is for a simulation, and it needs to be able to run potentially thousands to millions of times without falling over. The Dantzig solver I'm using right now, for instance, works in $O(n^3)$ as long as my $A$ matrix is positive semi-definite. If the $A$ matrix isn't, it's possible for it to diverge or take exponential time. Do I have similar guarantees with interior point methods?
-