eventsmodelicaopenmodelicasystemmodeler

How to use boolean switches in Modelica to prevent draining a stock below zero?

Problem description

I am building a library to support System Dynamics (SD) like modeling in Modelica. Unlike the freely available library by Cellier et al. I am convinced that one may well make use of acausal connectors: Transmitting the stock value as "potential" via connectors enables to build compact components (e.g. flows = processes).

In SD we would maybe distinguish material ("mass") stocks from informational stocks that can become negative. To support this I am using the following definitions for a mass port (here giving the definition for a StockPort - its the counterpar, a FlowPort, would simply have Boolean input variables instead of output variables and is given later on):

 connector StockPort "Used to represent stock and flow connections"
    Real stock "Current value of material in the stock";
    flow Real rate "Flow that affects the stock";
    // Boolean switches
    output Boolean stopInflow "True indicates that nothing can flow into the stock";
    output Boolean stopOutflow "True indicates that nothing can flow out of the stock";
 end StockPort;

The Boolean switches indicate for each port of a stock whether filling or draining is allowed.

For a "material stock" the stopOutflow switch should prevent the stock from being drained below zero. Unfortunately in the following example this will not work out: The stock will be drained just slightly below zero.

Minimal example using connectors

The following TestModel uses these building blocks:

The main model simply initiates a stock with initialValue = 5 that is connected to a process of linear decline with declineRate = 1.

model TestModel "Stop draining a stock below zero"

  function constrainedRate "Set rate for a port according to signals from stock" 
    input Real indicatedRate "Proposed rate for port of flow element";
    input Boolean stopInflow "Signal from connected stock";
    input Boolean stopOutflow "Signal from connected stock";
    output Real actualRate "The rate to use";
  protected
    // check whether indicated rate is negative (e.g. an inflow to the connected stock)
    Boolean indRateIsInflow = indicatedRate < 0;
  algorithm
    // set rate to zero if stopSignal matches character of flow
    actualRate := if indRateIsInflow and stopInflow 
          then 0 
        elseif not indRateIsInflow and stopOutflow 
          then 0 
        else indicatedRate;
  end constrainedRate;

  connector FlowPort "Used to represent stock and flow connections"
    Real stock "The current stock level (e.g. Potential) of a connected stock or flow data for special stocks";
    flow Real rate "Flows that affect the material stock";
    input Boolean stopInflow "True indicates that nothing can flow into the stock";
    input Boolean stopOutflow "True indicates that nothing can flow out of the stock";
  end FlowPort;

  connector StockPort "Used to represent stock and flow connections"
    Real stock "Current value of stock";
    flow Real rate "Flow that affects the stock";
    output Boolean stopInflow "True indicates that nothing can flow into the stock";
    output Boolean stopOutflow "True indicates that nothing can flow out of the stock";
  end StockPort;

  model MaterialStock "Stock that cannot be drained below zero"
    StockPort outflow;
    parameter Real initialValue;
  protected
    Real x(start = initialValue);
  equation
    // rate of change for the stock
    der(x) = outflow.rate;
    // ports shall have level information for stock
    outflow.stock = x;
    // inflow to stock is unrestricted
    outflow.stopInflow = false;
    // provide Boolean signal in case of negative stock
    outflow.stopOutflow = x <= 0;
  end MaterialStock;

  model LinearDecline "Decline of stock at a constant rate"
    FlowPort massPort;
    parameter Real declineRate(min = 0) "Rate of decline (positive rate diminishes stock)";
  protected
    // a positive rate should drain the stock (which here matches Modelica's rule for flows)
    Real rate(min = 0);
  equation
    rate = declineRate;
    // observe stock signals and constraints
    assert(rate >= 0, "Rate must be positive and will be set to zero", level = AssertionLevel.warning);
  // set the rate according to constraints given by stock
    massPort.rate = constrainedRate( max(rate, 0), massPort.stopInflow, massPort.stopOutflow );
  end LinearDecline;

  // main model
  MaterialStock stock( initialValue = 5 );
  LinearDecline process( declineRate = 1 );
equation
  connect( stock.outflow, process.massPort );
end TestModel;

Simulating the model using DASSL from StartTime = 0 to StopTime = 10 reveals the expected behavior for the variable stock.outflow.stock:

Stock Value

Unfortunately the value is slightly below zero at t = 5.0 and later.

Somehow the event (stock value <= 0) is detected too late. What can I do?

(So far the imo unelegant cure has been to use a when event (state-event) to reinit the stock value to zero. My experiments with using noEvent wrappers on if statements and Boolean conditions have also not been successful.)

Solution

  • Could you test this solution? Using Modelica.Constants.eps works for me. I also used Dassl and it worked for different stepsizes. I changed the following line from:

    // provide Boolean signal in case of negative stock
outflow.stopOutflow = x <= 0;

    to

    // provide Boolean signal in case of negative stock
outflow.stopOutflow = x <= Modelica.Constants.eps;

    Then, the output array(viewed in python for this post):

    Output using a zero instead of eps:
[ 5.00000000e+00  4.00000000e+00  3.00000000e+00  2.00000000e+00
  1.00000000e+00  0.00000000e+00 -7.37276906e-12 -7.37276906e-12
 -7.37276906e-12 -7.37276906e-12 -7.37276906e-12 -7.37276906e-12
 -7.37276906e-12 -7.37276906e-12]
Output using eps:
[5. 4. 3. 2. 1. 0. 0. 0. 0. 0. 0. 0. 0.]