4

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:

  • function constrainedRate( indicated rate, stopInflow, stopOutflow) is used to return a rate that will comply with the given constraints (i.e. boolean switches)

  • connector StockPort as described above

  • connector FlowPort the counterpart to the StockPort
  • model MaterialStock a stock-component with one StockPort that shall not be drainable below zero
  • model LinearDecline a flow-element with one FlowPort (i.e. a Sink) that models draining a connected stock at a constant rate (here set to 1)

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.)

gwr
  • 465
  • 6
  • 17
  • Have you tried setting the booleans not by comparing with 0 but with some small value like Modelica.Constants.eps or similiar? E.g.: stopOutflow = x <= Modelica.Constants.eps – f.wue Feb 18 '19 at 15:41
  • @f.wue Thanks. I just tried that, but unfortunately to no avail. – gwr Feb 18 '19 at 15:50

1 Answers1

2

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.]
f.wue
  • 837
  • 8
  • 15
  • I have tried DASSL (fixed interval sizes) and Euler (different step sizes) but your solution will not work out in either Wolfram System Modeler (4.3) or Open Modelica. I added `assert( x >= 0, "Stock must always be positive")` to have the simulation stop in case of negativity. It always stops. – gwr Feb 19 '19 at 09:29
  • Ok. I tested it in dymola, maybe that's the reason for the different results. Are the negative values on your machine also quite small (e.g. e-12) or bigger? – f.wue Feb 19 '19 at 09:33
  • What does work is having `outflow.stock = noEvent(if x > 0 then x else 0)`. In that case the numerical "issue" stays inside the component where if absolutely needed a `when x <= eps then reinit( x, 0 ); end when;` may take care of it. – gwr Feb 19 '19 at 09:42
  • Ok. In the end it depends what you want to model. If stock-values near zero occur often (with near zero i mean e-4 or so), you maybe need to use the noEvent or when solution. However, if you model bigger values, your model most likely is not valid for such small numbers. Therefore, you could use a bigger than eps (which is e-15 i believe). Something like <= 1e-5 should most definitly work with any solver/software in my opinion. If simulation time however is no constraint for you, use the noEvent method. – f.wue Feb 19 '19 at 10:23
  • I do not understand why modeling bigger values would invalidate the nonzero-checks - they simply would not apply? – gwr Feb 19 '19 at 10:27
  • What i am trying to say is: If your model is designed to apply for big values, you may also "cut-off" small values near zero. Sure, this is a simplification. However, compared to generally big numbers, cutting of 1e-5 does not provoke a big error. But your model gains efficiency as you don't rely on the reinit. If however try to model small values, cutting-off near zero will provoke greater erros and should therefore not be done. – f.wue Feb 19 '19 at 10:34
  • Thank you for the explanation. System Dynamics modeling is usally about social and economic systems so that small numbers would usually only appear as some fractions (e.g. 0.01 % would be maybe suffice for quite small). – gwr Feb 19 '19 at 10:42
  • Let us [continue this discussion in chat](https://chat.stackoverflow.com/rooms/188665/discussion-between-f-wue-and-gwr). – f.wue Feb 19 '19 at 10:50