Constrained Newton-Raphson root finding

Question

Original Question

I have a set of non-linear equations and I need to find the root where a subset of my solution vector is constrained to be greater than or equal to 0. I have implemented the Newton-Raphson algorithm but I am finding that some of the quantities that I need to remain positive are going negative. I am familiar with using Lagrange multipliers to enforce constraints in an optimization problem but I'm unsure about how to do the same for root finding.

It had occurred to me that I could map my variables which must be positive to log space (i.e., solve for $\log(x_i)$ rather than $x_i$) but since zero is an important case I was worried that this might cause numerical problems.

From a previous answer it looks like an interior point method might be appropriate but it would involve modifying my solver rather than just the residual equation which I would like to avoid.

Further Description

The equations are a series of non-linear equations arising from a material model which incorporates irrecoverable deformation ( $E^p$ ). The way that this is typically achieved is through internal state variables ($\xi$) which are governed through some evolution equations which are determined by some evolution rate ( $\dot{\gamma}$ ).

The evolution equations take the form of:

$\dot{E^p} = \dot{E}\left(E^p,\xi, \dot{\gamma}\right)$

$\dot{\xi} = \dot{\xi}\left( \xi, \dot{\gamma}\right)$

Subject to some additional onset condition ( yield function ):

$F=F(E - E^p)$

where $E$ is the total deformation which is a function which is <= 0. The problem is subject to a Kuhn-Tucker condition which can be summarized as:

$\dot{\gamma}F = 0$

Implying that the evolution rate is zero if $F < 0$ and $\dot{\gamma}>0$ if $F=0$.

My unknown vector is something like:

$x=\left[ E^p, \xi, \gamma \right]$

My residual is something like:

$R=\left[{ E }^{p,expected} - E^p, \xi^{expected} - \xi, \dot{\gamma}F + \langle F \rangle \right]$

Where the values marked $\left(\cdot\right)^{expected}$ are the results of the evolution equations and $\langle \cdot \rangle$ are the Macaulay brackets which are defined as:

$\langle x \rangle = \frac{1}{2}\left( x + abs( x ) \right)$

I include the Macaulay term because, otherwise, the case of F > 0 and $\dot{\gamma}=0$ would be accepted which is incorrect. I tried using if statement in the residual calculation but that caused other issues and this seems to be performing better.

Answer 1

I thought I would come back to this now that I have an answer that works. The problem was multifaceted and I tried several solutions that worked to varying degrees before eventually settling on a solution.

Things I tried that didn't work on their own:

1. Homotopy solver that tried to implement barriers

   $g( x, s ) = (1 - \frac{ 1 } { a( s ) } ) R( x ) + \frac{ 1 }{ a( s ) } b( x, s )$

   where

   $b( x, s ) = e^{ s a ( b - x ) } - 1$

   $a( s ) = e^{ A s }$

   Where $R$ is the original residual, $s$, is the pseudo time, and $A$ is some value large enough to make sure the barrier holds. I found that this gave a pretty smooth variation from the barrier function to the true function if A is at least 10 or so. This worked for toy problems but not for my actual problem.

2. Newton homotopy solver:

   $g( x, s ) = R( x ) + ( 1 - s ) R( x_0 )$

   I like this homotopy and ended up using it for my final non-linear equation solve. In the solve I first try s = 1 and then cutback if required.

3. Performing multiple nested Newton-Raphson solves. Instead of trying to solve for $E^p$, $\xi$, and $\dot{\gamma}$ simultaneously I solved for $\xi$ analytically and $E^p$ numerically for a given value of $\dot{\gamma}$ as each iteration to solve for $\dot{\gamma}$. This dramatically helped with my convergence for many problems but didn't solve the more difficult cases.

4. Worked through my code and found a bug which was likely causing unintended behavior. Again, this helped, but didn't solve my problems.

5. Tried limiting the line-search in the Newton Raphson to always be in the correct domain. This helped, but sometimes didn't allow for convergence.

Now for what actually worked:

I had unintentionally set up my residual equation such that if the error was too great it would pass over into a new basin of attraction which pulled it in the wrong direction. I needed to modify my residual equation to prevent this behavior. Note that in the original residual for $\dot{\gamma}$

$R^{\dot{\gamma}} = \langle F \rangle + \dot{ \gamma } F$

The problem with this is that if $F$ gets very large, eventually the Jacobian can go positive due to the second term which will try to drive $\dot{\gamma}$ negative. So in the end I needed to:

1. Use a Newton homotopy solver ( this isn't necessarily required right now but helps when the problem gets stiff ).
2. Break the monolithic solve into sub-problems.
3. Modify my residual to exclude the problematic behavior.

In the context of the original equation, I believe 3 is probably the crucial point.