Logo Search packages:      
Sourcecode: gnucap version File versions  Download package

inverse.cc

/*$Id: inverse.cc,v 24.5 2003/04/27 01:05:05 al Exp $ -*- C++ -*-
 * Copyright (C) 2002 Telford Tendys
 * Author: Telford Tendys <telford@triode.net.au>
 *
 * This file is part of "Gnucap", the Gnu Circuit Analysis Package
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 2, or (at your option)
 * any later version.
 *
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA
 * 02111-1307, USA.
 *------------------------------------------------------------------
 * 
 * The purpose is to find the inverse of a function by numerical
 * methods. We are given:
 *      f(x) = y (where y is constant and known)
 *      f(x) may be evaluated as a black-box for any x
 *      f'(x) is also evaluated for any x
 *
 * We need to find x (starting from initial guess x0)
 * 
 */
#define NDEBUG

#include <math.h>
#include <assert.h>
#include "e_elemnt.h"
#include "u_opt.h"

#ifdef TEST_MODE
#include <stdio.h>
#endif // TEST_MODE

/*
 * FIXME: Some arbitrary epsilon constants have been used here
 */
static void newtons_method_of_guessing( ELEMENT *E, double delta_y )
{
#ifdef TEST_MODE
      printf( "Using Newton's Method\n" );
#endif // TEST_MODE
      if( fabs( E->_y0.f1 ) > 1e-50 )
      {
            E->_y0.x -= delta_y / E->_y0.f1;
      }
      else
      {
            if( E->_y0.f1 > 0 )
                  E->_y0.x -= delta_y * 1e50; /* This is not good */
            else
                  E->_y0.x += delta_y * 1e50; /* This is not good either */
      }
}

/*
 * Inputs:
 *    y is in E->_y0.f0
 *    x0 (initial guess) is in E->_y0.x
 *
 * E->tr_eval() will evaluate new values for E->_y0.f0 and E->_y0.f1
 *
 * Outputs:
 *    y goes back into E->_y0.f0 (actually our best result goes here)
 *    x goes into E->_y0.x
 *
 * Return values:
 *    0 indicates successful convergence
 *   -1 indicates ran out of iterations
 *   -2 indicates problem with E
 *
 * (probably should have some NaN checks that it can bomb out with)
 */
int ELEMENT::inverse()
{
      int iter, state = 0;
      double target_y, scale;
      double low_x = 0, low_y = 0;
      double high_x = 0, high_y = 0;
      register double delta_y; /* A workhorse variable */

      if( !has_tr_eval()) return( -2 ); /* No function available */
      target_y = _y0.f0;
      scale = fabs( target_y );
      if( scale < 1e-50 ) scale = 1e-50; /* FIXME: use options */
      for( iter = 200; iter--; )
      {
#ifdef TEST_MODE
            printf( "Before eval: state=%d target_y=%g x=%g f0=%g f1=%g\n",
                        state, target_y, _y0.x, _y0.f0, _y0.f1 );
#endif // TEST_MODE
            common()->tr_eval(this);
            delta_y = _y0.f0 - target_y;
#ifdef TEST_MODE
            printf( "After eval: delta_y=%g x=%g f0=%g f1=%g\n",
                        delta_y, _y0.x, _y0.f0, _y0.f1 );
#endif // TEST_MODE
            if( fabs( delta_y ) / scale < 1e-10 ) return( 0 ); /* FOUND! */
            switch( state )
            {
/*
 * Nothing is valid, we have an initial guess only
 * Figure out what side of the fence it sits on and thus get
 * either a valid high or a valid low. From there, use a 
 * Newton's Method step to have a stab at a second guess.
 * Frequently this stab will actually be very good but we
 * don't actually trust it.
 */
                  case 0:
                        low_x = high_x = _y0.x;
                        low_y = high_y = _y0.f0;
                        if( delta_y > 0 )
                        {
                              state = 1;
                        }
                        else
                        {
                              state = 2;
                        }
                        newtons_method_of_guessing( this, delta_y );
                        break;
/*
 * If possible, straddle a zero crossing and jump to binary search
 * but if we can't find the straddle then use Newton's Method as a fallback.
 */
                  case 1:
                        assert( high_y > target_y );
                        assert( low_x == high_x );
                        assert( low_y == high_y );
                        if( delta_y < 0 )  /* Is he straddlin'? */
                        {
                              if( high_x > _y0.x )
                              {
                                  low_x = _y0.x;
                                    low_y = _y0.f0;
                                    state = 3;
                              }
                              else
                              {
                                  high_x = _y0.x;
                                    high_y = _y0.f0;
                                    state = 4;
                              }
 midpoint_binary_search:
#ifdef TEST_MODE
                              printf( "Using Binary Search\n" );
#endif // TEST_MODE
                              _y0.x = ( low_x + high_x ) / 2.0;
                              break;
                        }
                        if( high_x > _y0.x )
                        {
                              low_y = _y0.f0;
                              low_x = _y0.x;
                        }
                        else
                        {
                              high_y = _y0.f0;
                              high_x = _y0.x;
                        }
                        newtons_method_of_guessing( this, delta_y );
                        state = 5;
                        break;
/*
 * If possible, straddle a zero crossing and jump to binary search
 * but if we can't find the straddle then use Newton's Method as a fallback.
 */
                  case 2:
                        assert( low_y < target_y );
                        assert( low_x == high_x );
                        assert( low_y == high_y );
                        if( delta_y > 0 )  /* Is he straddlin'? */
                        {
                              if( low_x < _y0.x )
                              {
                                  high_x = _y0.x;
                                    high_y = _y0.f0;
                                    state = 3;
                              }
                              else
                              {
                                  low_x = _y0.x;
                                    low_y = _y0.f0;
                                    state = 4;
                              }
                              goto midpoint_binary_search;
                        }
                        if( high_x > _y0.x )
                        {
                              low_y = _y0.f0;
                              low_x = _y0.x;
                        }
                        else
                        {
                              high_y = _y0.f0;
                              high_x = _y0.x;
                        }
                        newtons_method_of_guessing( this, delta_y );
                        state = 6;
                        break;
/*
 * We have a straddle, narrow the region to a smaller straddle
 * and attempt to apply Newton's Method...
 * if that fails use a binary search as a fallback
 */
                  case 3:
                        assert( high_y > target_y );
                        assert( low_y < target_y );
                        assert( high_x > low_x );
                        assert( _y0.x >= low_x );
                        assert( _y0.x <= high_x );

                        if( delta_y > 0 )
                        {
                              high_x = _y0.x;
                              high_y = _y0.f0;
                        }
                        else
                        {
                              low_x = _y0.x;
                              low_y = _y0.f0;
                        }
                        newtons_method_of_guessing( this, delta_y );
                        if( _y0.x < low_x ) goto midpoint_binary_search;
                        if( _y0.x > high_x ) goto midpoint_binary_search;
                        break;
/*
 * We have a straddle, narrow the region to a smaller straddle
 * and attempt to apply Newton's Method...
 * if that fails use a binary search as a fallback
 */
                  case 4:
                        assert( high_y < target_y );
                        assert( low_y > target_y );
                        assert( high_x > low_x );
                        assert( _y0.x >= low_x );
                        assert( _y0.x <= high_x );

                        if( delta_y < 0 )
                        {
                              high_x = _y0.x;
                              high_y = _y0.f0;
                        }
                        else
                        {
                              low_x = _y0.x;
                              low_y = _y0.f0;
                        }
                        newtons_method_of_guessing( this, delta_y );
                        if( _y0.x < low_x ) goto midpoint_binary_search;
                        if( _y0.x > high_x ) goto midpoint_binary_search;
                        break;
/*
 * We still haven't found a straddle and we have already had two tries
 * at Newton's Method. It we have a straddle by now then fine,
 * if not then start to expand our span.
 */
                  case 5:
                        assert( high_y > target_y );
                        assert( low_y > target_y );
                        assert( high_x > low_x );
                        if( delta_y < 0 )  /* Is he straddlin'? */
                        {
                              if( high_x > _y0.x )
                              {
                                  low_x = _y0.x;
                                    low_y = _y0.f0;
                                    state = 3;
                              }
                              else
                              {
                                  low_x = high_x;
                                    low_y = high_y;
                                  high_x = _y0.x;
                                    high_y = _y0.f0;
                                    state = 4;
                              }
                              goto midpoint_binary_search;
                        }
#ifdef TEST_MODE
                        printf( "Widening the span in search of straddle\n" );
#endif // TEST_MODE
                        if( _y0.x < low_x )
                        {
                              low_y = _y0.f0;
                              low_x = _y0.x;
                              _y0.x = high_x + high_x - low_x;
                        }
                        else if( _y0.x > high_x )
                        {
                              high_y = _y0.f0;
                              high_x = _y0.x;
                              _y0.x = low_x + low_x - high_x;
                        }
                        else
                        {
                              _y0.x = low_x + low_x - high_x;
                        }
                        break;
/*
 * Still no straddle... same as case 5 except we are on the other side
 */
                  case 6:
                        assert( high_y < target_y );
                        assert( low_y < target_y );
                        assert( high_x > low_x );
                        if( delta_y > 0 )  /* Is he straddlin'? */
                        {
                              if( low_x < _y0.x )
                              {
                                  high_x = _y0.x;
                                    high_y = _y0.f0;
                                    state = 3;
                              }
                              else
                              {
                                  low_x = _y0.x;
                                    low_y = _y0.f0;
                                    state = 4;
                              }
                              goto midpoint_binary_search;
                        }
#ifdef TEST_MODE
                        printf( "Widening the span in search of straddle\n" );
#endif // TEST_MODE
                        if( _y0.x < low_x )
                        {
                              low_y = _y0.f0;
                              low_x = _y0.x;
                              _y0.x = high_x + high_x - low_x;
                        }
                        else if( _y0.x > high_x )
                        {
                              high_y = _y0.f0;
                              high_x = _y0.x;
                              _y0.x = low_x + low_x - high_x;
                        }
                        else
                        {
                              _y0.x = high_x + high_x - low_x;
                        }
                        break;
            }
      }
/*
 * We ran out of iterations, there isn't much we can do here
 * so just return an error. Possibly the error code should be
 * based on the state but at the moment the error code doesn't
 * give much info about how close we got
 */
      return( -1 );
}

/*
 * Telford's Knee-chord method of estimating a linear model for saturated inductors
 *
 * NB: All this is written from the point of view of an inductor (see d_coil.cc)
 * but the principle is not limited to any particular device.
 *
 * -- Find the point that is off the main curve which is the system operating point
 *    outside of this device (if we presume that we are the only non-linear device
 *    in the world then this point is the intersection of the load-line of everything
 *    except this device and the previous estimated linear model of this device).
 *
 *    The point in question (call it OP1) is:
 *          (current=_y0.x, flux={_y1.f0+_y1.f1*(y0.x-y1.x)})
 *
 *    This is obtained by taking the most recent current estimate and plugging it into
 *    the PREVIOUS linear estimate of the device behaviour.
 *
 * -- Draw a line horizontal from OP1 until it touches the non-linear curve
 *    for this device (this is standard inversion of a well-behaved function,
 *    we use a hybrid of Newton's Method and a binary search for this purpose).
 *    By horizontal, I mean keep constant flux and adjust the current.
 *
 *    This leads us to OP2 which is:
 *          (current=<result of function inversion>, flux=<same as OP1>)
 *
 * -- Draw a vertical line from OP1 to touch the non-linear curve
 *    (this is just a standard tr_eval() which has just been done so we merely
 *    save the _y0.f0 and _y0.x values and restore them after we have finished
 *    with function inversions). Call this point OP3 (this is the operating point
 *    that the standard Newton's Method uses as a working estimate of the system,
 *    the only difference between the Knee-chord method and the standard method is
 *    the value of _y0.f1)
 *
 * -- Find the slope of the interval connecting OP2 and OP3 and use this slope
 *    as the _y0.f1 value rather than the simple derivative of the curve.
 *
 *    Using this slope has a number of advantages over Newton's Method with regard
 *    to stability. If OP1 is close to the non-linear curve then we are close
 *    to convergence and we can say that the non-linear curve is close to linear
 *    in the local region. Thus OP2 and OP3 will be close together and the inverval
 *    between them will essentially be the same as the local function derivative
 *    (actually, once OP2 and OP3 are sufficiently close, it is better to switch
 *    back to Newton's Method since the derivative will be calculated more accurately).
 *    If OP2 and OP3 are far apart then the slope becomes a broad average that
 *    estimates the general behaviour of the non-linear curve in the region near
 *    where we are operating. This is important when we are far from converged
 *    so that stability is maintained.
 *
 * Probably this should return a success/failure code but I don't know what
 * to use the returned value for so don't bother right now
 */
void ELEMENT::knee_chord()
{
      double save_x, save_f0, save_f1;

      if( OPT::strategy != stKNEECHORD ) return; /* This strategy not in use */
      if( !has_tr_eval()) return; /* No function available */
/*
 * Save the standard values -- these are OP3
 */
      save_x = _y0.x;
      save_f0 = _y0.f0;
/*
 * Save this too, just in case we can't get the result we want
 * (fallback to standard algorithm)
 */
      save_f1 = _y0.f1;
/*
 * Project the old linear estimate of device behaviour onto the new coil current
 * in order to come up with a flux linkage value (in _y0.f0) gives us OP1
 */
      _y0.f0 = _y1.f0 + _y1.f1 * ( _y0.x - _y1.x );
/*
 * Invert the function so that we find a _y0.x value that gives _y0.f0
 * (thus we get OP2)
 */
      if( 0 > inverse())
      {
            untested();
            error( bPICKY, long_label() + ": function inversion failed\n" );
            _y0.f1 = save_f1; /* Fallback */
      }
      else
      {
/*
 * Find the all-important chord gradient
 */
            double numerator = _y0.f0 - save_f0;
            double denominator = _y0.x - save_x;

            if( fabs( denominator ) > 1e-50 )
            {
                  _y0.f1 = numerator / denominator;
#ifdef TEST_MODE
                  printf( "Chord gradient is %g\n", _y0.f1 );
#endif // TEST_MODE
            }
            else
            {
#ifdef TEST_MODE
                  printf( "Fallback because OP2 and OP3 are too close\n" );
#endif // TEST_MODE
                  _y0.f1 = save_f1; /* Fallback */
            }
      }
/*
 * Back to the original point again
 * (Come to think of it, leaving it at OP2 might also work)
 */
      _y0.x = save_x;
      _y0.f0 = save_f0;
}

Generated by  Doxygen 1.6.0   Back to index