[SciPy-User] SCIPY FSOLVE

[josef.pktd at gmail.com]

On Mon, Oct 12, 2009 at 7:55 AM, Etrade Griffiths
<etrade.griffiths at dsl.pipex.com> wrote:
> Hi
>
> I am trying to solve directly a series of equations describing flow
> in a network using FSOLVE but have not had much success so far.  A
> small example is given below.
>
> I have a feeling that the problem may be ill-conditioned because
> there are two sources of small numbers in the equations: one
> associated with small coefficients and the other associated with
> small differences in pressure between adjacent nodes.  The main
> problem seems to be that FSOLVE ends up with estimates of the
> variables that result in raising a negative number to a power, so
> that some of the residuals become -1.#IND in value.  I tried trapping
> this type of error but FSOLVE does not appear to be able to move
> towards the solution.  I would be grateful for any suggestions on how
> to "encourage" FSOLVE.
>
> Thanks in advance
>
> Alun Griffiths
>
> # =========================
> #
> # Simple model of gathering network
> # Uses SCIPY FSOLVE to solve for network pressures directly
> #
>
> import math
> import scipy.optimize
>
> # Set up the system of network equations
>
> def NetworkEquations(x):
>
>     # Transfer current guess to variables
>
>     p2 = x[0]
>     q1 = x[1]
>     q2 = x[2]
>     q3 = x[3]
>
>     # Define constants
>
>     c = [300.0, 2.00E-5, 1.50E-5]
>     n = [0.50, 0.75, 0.70]
>
>     p1 = 200.0
>     p3 = 2000.0
>
>     # Define the residuals
>
>     residuals = []
>
>     # ... Kirchoff
>
>     residuals.append( q1 - q2 - q3 )
>
>     # ... pressure drop along pipeline
>
>     curr_err = q1 - (  c[0] * ( p2 ** 2 - p1 ** 2 ) ) ** n[0]
>     residuals.append(curr_err)
>
>     # ... pressure drop over well 1
>
>     curr_err = q2 - c[1] * ( p3 ** 2 - p2 ** 2 ) ** n[1]
>     residuals.append(curr_err)
>
>     # ... pressure drop over well 2
>
>     curr_err = q3 - c[2] * ( p3 ** 2 - p2 ** 2 ) ** n[2]
>     residuals.append(curr_err)
>
>     # All equations evaluated so return residuals
>
>     return residuals
>
>
> # Solve equations using Broyden's method
>
> x0 = [250.0, 2.0, 1.0, 1.0]
> x1 = scipy.optimize.fsolve(NetworkEquations,x0)
> print x1
>
>
> _______________________________________________
> SciPy-User mailing list
> SciPy-User at scipy.org
> http://mail.scipy.org/mailman/listinfo/scipy-user
>

I never used optimize.fsolve, but in your case, I would try
to reparameterize p2 so it has to be between between
p1 and p3,

p1<=p2<=p3

maybe as a fraction in interval [0,1] and then transform it
to the interval [p1,p3]

There might be other ways to impose the constraint, but
my first attempt is usually reparameterization.

Josef