Package ptolemy.domains.algebraic.kernel

Class AlgebraicLoopDirector

  • All Implemented Interfaces:
    java.lang.Cloneable, Executable, Initializable, Changeable, Debuggable, DebugListener, Derivable, ModelErrorHandler, MoMLExportable, Moveable, Nameable
    public class AlgebraicLoopDirector
extends StaticSchedulingDirector
    A director for solving algebraic loops. This director initializes all inputs that have a defaultValue parameter, and then executes the model according to the specified solver method until all inputs change by less than a specified errorTolerance. The methods implemented are:
    1. SuccessiveSubstitution: This simple strategy executes the model until all inputs have converge.
    2. NewtonRaphson: The Newton-Raphson iteration.
    3. Homotopy: A homotopy method that uses Newton-Raphson iteration and Euler updates that search along a curve on which a linear approximation to the problem is successively replaced with the actual function.
    In all cases, the number of iterations is limited to maxIterations, and an exception will be thrown if convergence hasn't occurred by then.

    The errorTolerance may be given as a parameter to an individual port (just add a parameter named "errorTolerance" to the port). For any port for which there is no such parameter, the errorTolerance parameter of this director will be used. Note that if the errorTolerance of a port is changed during a run, the new value is ignored until the next run.

    In all cases, the problem being solved has the form: 

         -----------
     |  -----  |
     -->| g |---
     x  -----
    where x is initially the vector of values corresponding to input ports that have a non-null defaultValue parameter, and g is the network of actors connecting these ports.

    This class solves an algebraic loop of the form x=g(x) 

         -----------
     |  -----  |
     -->| g |---
     x  -----
    where x is initially the vector of values corresponding to input ports that have a non-null defaultValue parameter, and g is the network of actors connecting these ports.

    For the SuccessiveSubstitution, we simply evaluate g repeatedly until either all signals do not change by more than the errorTolerance, or until there have been maxIterations. Note that it is possible to reach a fixed point where some or all signals are infinite.

    For NewtonRaphson, we solve for g(x)=x by solving f(x)=0, where f(x) = x-g(x). This is done by iterating as follows: 

       x_n+1 = x_n - f(x_n)/f'(x_n) = x_n - (x_n - g(x_n))/(1 - g'(x_n)) .
    To estimate g'(x_n), we do 
       g'(x_n) = (g(x_n + d) - g(x_n))/d
    where d is the delta parameter.

    For Homotopy, we solve f(x) = x - g(x). The problem is reformulated as H(s, lambda, x0) = s - lambda (g(s+x0)-x0), where lambda has an initial value of 0 and is successively increased to 1, s is a coordinate transformation defined so that x = s+x0, where x0 is the initial iterate.
    The implementation is equal to Program 3 of Eugene L. Allgower and Kurt Georg, Introduction to Numerical Continuation Methods, Classics in Applied Mathematics, Vol. 45, SIAM, 2003. However, the implementation by Allgower and Georg assumes an initial iterate of 0, which is the reason for the above coordinate transformation.

    FIXME: Questions:

    • This implementation checks for convergence of all input ports, not just those with defaultValue parameters (i.e., not just x). Is this what we want?
    • This implementation uses the defaultValue on every firing. Should the default value instead be replaced by the previous solution on all but the first iteration? FIXME: Yes, this will need to be done as it increases robustness and efficiency.
    • Should delta be able to be different for each variable? I.e., should each input port that has a defaultValue parameter also be able to have a delta parameter? FIXME: Yes, the magnitudes of variables can be a few orders of magnitude different, e.g., if one port as a control signal, and another an air pressure in Pascal.
    • Instead of estimating g'(x_n), we want to be able to optionally identify g'(x_n). But this is really tricky when x_n is a vector (i.e. where there are multiple input ports with defaultValue). FIXME: This is addressed in the Homotopy algorithm which uses Broyden updates to estimate the Jacobian. A similar mechanism could be added to the Newton algorithm in future work.
    • This code is not at all optimized and may be quite inefficient.
    Since:
    Ptolemy II 10.0
    Version:
    $Id$
    Author:
    Edward A. Lee and Michael Wetter
    Pt.AcceptedRating:
    Red (pwhitake)
    Pt.ProposedRating:
    Yellow (eal)

    • Field Detail

      • errorTolerance

        public Parameter errorTolerance
        The default tolerance for determining when convergence has occurred. When the current value of an input port differs from the previous value by less than the errorTolerance, then we declare it to hfave converged. This parameter gives a default value that will be used if the port has no parameter named "errorTolerance". This is a double with default value 1E-4.

      • maxIterations

        public Parameter maxIterations
        The maximum number of allowed iterations before the director will declare a failure to converge. This is an integer that defaults to 1000.

      • method

        public StringParameter method
        The method to be used to solve algebraic loops. This is a string that is one of "Homotopy" (the default), "NewtonRaphson" or "SuccessiveSubstitution".

      • _breakVariables

        protected java.util.List<AlgebraicLoopReceiver> _breakVariables
        The list of receivers for all break variables.

      • _g_n

        protected double[] _g_n
        Current value of the loop function g(x_n).

      • _nVars

        protected int _nVars
        Number of break variables.

      • _tolerance

        protected double[] _tolerance
        Tolerance for each iteration variable.

      • _x_n

        protected double[] _x_n
        Current value of the iteration variables x_n.

    • Constructor Detail

      • AlgebraicLoopDirector

        public AlgebraicLoopDirector​(CompositeEntity container,
                             java.lang.String name)
                      throws IllegalActionException,
                             NameDuplicationException
        Construct a director in the given container with the given name. The container argument must not be null, or a NullPointerException will be thrown. If the name argument is null, then the name is set to the empty string. Increment the version number of the workspace.
        Parameters:
        container - Container of the director.
        name - Name of this director.
        Throws:
        IllegalActionException - If the director is not compatible with the specified container.
        NameDuplicationException - If the name collides with an attribute in the container.

    • Method Detail

      • fire

        public void fire()
          throws IllegalActionException
        Prefire and fire actors in the order given by the scheduler until the iteration converges. An iteration converges when a pass through the schedule does not change the status of any receiver.
        Specified by:
        fire in interface Executable
        Overrides:
        fire in class StaticSchedulingDirector
        Throws:
        IllegalActionException - If an actor violates the monotonicity constraints, or the prefire() or fire() method of the actor throws it.

      • _evaluateLoopFunction

        protected void _evaluateLoopFunction​(double[] x,
                                     double[] g)
                              throws IllegalActionException
        Evaluate the loop function for x and save the result in g. This function is called by the solver to evaluate the loop function.
        Parameters:
        x - Input to the loop function.
        g - Double vector of the same size as x. The result will be stored in this function.
        Throws:
        IllegalActionException - If the prefire() method returns false having previously returned true in the same iteration, or if the prefire() or fire() method of the actor throws it, or if evaluating the function yields a value that is not a double.

      • newReceiver

        public Receiver newReceiver()
        Return a new FixedPointReceiver. If a subclass overrides this method, the receiver it creates must be a subclass of FixedPointReceiver, and it must add the receiver to the _receivers list (a protected member of this class).
        Overrides:
        newReceiver in class Director
        Returns:
        A new FixedPointReceiver.

      • _clearAllDestinationReceivers

        protected void _clearAllDestinationReceivers​(Actor actor)
                                      throws IllegalActionException
        Call the send(index, null) method of each output port of the specified actor.
        Parameters:
        actor - The actor.
        Throws:
        IllegalActionException - If thrown while getting the width of a port, determining if a port is known or while sending data.

      • _fireActor

        protected void _fireActor​(Actor actor)
                   throws IllegalActionException
        Fire an actor. Call its prefire() method, and if that returns true, call its fire() method.
        Parameters:
        actor - The actor to be fired.
        Throws:
        IllegalActionException - If the prefire() method returns false having previously returned true in the same iteration, or if the prefire() or fire() method of the actor throws it.