Refactoring stanify: Mismatching initialization stock value

[hyunjimoon]

Stanify(.mdl) writes the following as production_start_rate in ode_vensim_function function in function block.

max(0, 1 + 2 + 2 * 1 - 2 + 2 * 1 / 7) 
+ 6 * max(0, 1 + 2 + 2 * 1 - 2 + 2 * 1 / 7) 
- 6 * max(0, 1 + 2 + 2 * 1 - 2 + 2 * 1 / 7) / 3

When wip_adjustement_time = 3, manufacturing_cycle_time = 6, inventory_adjustment_time = 8, time_to_average_order_rate = 8:

max(0, 1 + 2 + 2 * 1 - 2 + 2 * 1 / 8) 
+ 6 * max(0, 1 + 2 + 2 * 1 - 2 + 2 * 1 / 8) 
- 6 * max(0, 1 + 2 + 2 * 1 - 2 + 2 * 1 / 8) / 3;

When wip_adjustement_time = 2, manufacturing_cycle_time = 8, inventory_adjustment_time = 8, time_to_average_order_rate = 8:

max(0,
+ 8 * max(0, 1 + 2 + 2 * 1 - 2 + 2 * 1 / 8) 
- 8 * max(0, 1 + 2 + 2 * 1 - 2 + 2 * 1 / 8) 
      + max(0, 1 + 2 + 2 * 1 - 2 + 2 * 1  / 8) / 2)

Whereas manual computing replacing the following equation with values assigned in vensim (below) returns:

max(0, 
        (+ 8 * max(0, 1 + (+ (2 + 2) * 2 
                           - (2 + 2) * 2
                           )/8) 
         - 8 * max(0, 1 + (+ (2 + 2) * 2 
                           - (2 + 2) * 2
                           )/8) 
        )/8 
+ 
max(0, 
              1 + (+ 1 * (2 + 2) 
                   - 1 * (2 + 2)
                )/8
            )
)

production_start_rate_stock_init is computed using topologically-sorted abstract syntax tree: INIT: production_start_rate_stock_init = production_start_rate EQN:

• desired_production_start_rate = adj_for_wip + desired_production
• adj_for_wip= (desired_wip - work_in_process_inventory)/wip_adjustment_time

  
  = max(0, desired_production_start_rate) 
  = max(0, 
             adj_for_wip 
             + desired_production
    ) 

= max(0, (desired_wip - work_in_process_inventory) /wip_adjustment_time

• max(0, customer_order_rate + production_adjustment_from_inventory) )

  INIT: `work_in_process_inventory_init` = `desired_wip`
  EQN:` production_adjustment_from_inventory` = (desired_inventory - inventory)/inventory_adjustment_time

  = max(0, (desired_wip - desired_wip) /wip_adjustment_time

  • max(0, customer_order_rate + (desired_inventory - inventory)/inventory_adjustment_time ) )

    INIT: `inventory_init` = `desired_inventory`
    EQN: 
    `desired_wip` = `manufacturing_cycle_time` * `desired_production` 
    `desired_production` = `MAX(0, expected_order_rate + production_adjustment_from_inventory)`

    = max(0, (+ manufacturing_cycle_time * desired_production

    • manufacturing_cycle_time * desired_production )/wip_adjustment_time
      • max(0, customer_order_rate + (desired_inventory - desired_inventory )/inventory_adjustment_time ) )

        EQN: 
        `desired_inventory` = `MAX(0, desired_inventory_coverage + expected_order_rate)`

        = max(0, (+ manufacturing_cycle_time * desired_production

    • manufacturing_cycle_time * desired_production )/wip_adjustment_time
      • max(0, customer_order_rate + (max(0, desired_inventory_coverage + expected_order_rate)
      • max(0, desired_inventory_coverage + expected_order_rate) )/inventory_adjustment_time ) )

        
        INIT
        `expected_order_rate_init` = customer_order_rate

EQN: desired_production = MAX(0, expected_order_rate + production_adjustment_from_inventory) desired_wip_init = manufacturing_cycle_time * desired_production desired_inventory_init = expected_order_rate * desired_inventory_coverage

    = max(0, 
        (+ manufacturing_cycle_time * max(0, expected_order_rate + production_adjustment_for_inventory) 
          - manufacturing_cycle_time * max(0, expected_order_rate + production_adjustment_from_inventory)
        )/wip_adjustment_time
        + max(0, 
            customer_order_rate + 
           (max(0, `desired_inventory_coverage`+ customer_order_rate) 
           -  max(0, `desired_inventory_coverage` + customer_order_rate)
                   )/inventory_adjustment_time
           )
      )

INIT again from newly created (customer_order_rate always the second term lagging behind the first): expected_order_rate_init = customer_order_rate EQN:

• desired_inventory_coverage = minimum_order_processing_time + safety stock coverage

  = max(0, 
      (+ manufacturing_cycle_time * max(0, expected_order_rate + (desired_inventory - inventory)/inventory_adjustment_time) 
        - manufacturing_cycle_time * max(0, expected_order_rate + (desired_inventory - inventory)/inventory_adjustment_time)
      )/wip_adjustment_time
      + max(0, 
            customer_order_rate + 
             (max(0, + expected_order_rate * (minimum_order_processing_time + safety_stock_coverage) 
                 - max(0, expected_order_rate * (minimum_order_processing_time + safety_stock_coverage)
              )/inventory_adjustment_time
          )
    )

  INIT: expected_order_rate_init = customer_order_rate EQUATION: customer_order_rate = desired_inventory_coverage * expected_order_rate

  = max(0, 
      (+ manufacturing_cycle_time * max(0, customer_order_rate + (desired_inventory - desired_inventory)/inventory_adjustment_time) 
       - manufacturing_cycle_time * max(0, customer_order_rate + (desired_inventory - desired_inventory)/inventory_adjustment_time)
      )/wip_adjustment_time
      + max(0, 
            customer_order_rate + 
             (max(0,(+ customer_order_rate * (minimum_order_processing_time + safety_stock_coverage) 
                 -  max(0, customer_order_rate * (minimum_order_processing_time + safety_stock_coverage)
              )/inventory_adjustment_time
          )
    )

The following is numeric tracking using the values

• customer_order_rate = 1
• wip_adjustement_time = 3
• minimum_order_processing_time = 5
• manufacturing_cycle_time = 6
• inventory_adjustment_time = 7
• time_to_average_order_rate = 8
• safety_stock_coverage = 2

INSERT INIT production_start_rate_stock_init = production_start_rate

= max(0, desired_production_start_rate) 
= max(0, 
               adj_for_wip 
               + desired_production
      ) 

= max(0, 
             (desired_wip - work_in_process_inventory)
             /wip_adjustement_time
            + max(0, customer_order_rate + production_adjustment_from_inventory) 
      )
INSERT INIT work_in_process_inventory_init desired_wip

    = max(0, 
                 (desired_wip - desired_wip)
                 /3
                 + max(0, 1 + (desired_inventory - inventory)/inventory_adjustment_time
            ) 
      )
INSERT INIT inventory_init = desired_inventory

    = max(0, 
        (+ manufacturing_cycle_time * desired_production 
         - manufacturing_cycle_time * desired_production
        )/3
        + max(0, 
            1 + (desired_inventory - desired_inventory
                   )/7
           )
       )
INSERT INIT desired_wip_init = manufacturing_cycle_time * desired_production ,desired_inventory_init = expected_order_rate * desired_inventory_coverage

    = max(0, 
        (+ 6 * max(0, expected_order_rate + adjustment_for_inventory) 
          - 6 * max(0, expected_order_rate + adjustment_for_inventory)
        )/3
        + max(0, 
              1 + (+ expected_order_rate * desired_inventory_coverage 
                      - expected_order_rate * desired_inventory_coverage
                  )/7
            )
      )
INSERT expected_order_rate_init = customer_order_rate = 1

    = max(0, 
        (+ 6 * max(0, 1 + (desired_inventory - inventory)/inventory_adjustment_time) 
         - 6 * max(0, 1 + (desired_inventory - inventory)/inventory_adjustment_time)
        )/3
        + max(0, 
              1 + (+ 1 * (minimum_order_processing_time + safety_stock_coverage) 
                   - 1 * (minimum_order_processing_time + safety_stock_coverage)
                )/7
            )
      )
INSERT INIT

    = max(0, 
        (+ 6 * max(0, 1 + (desired_inventory - desired_inventory)/7) 
         - 6 * max(0, 1 + (desired_inventory - desired_inventory)/7)
        )/3
        + max(0, 
              1 + (+ 1 * (2 + 2) 
                   - 1 * (2 + 2)
                )/7
            )
      )

    = max(0, 
        (+ 6 * max(0, 1 + (+ desired_inventory_coverage * expected_order_rate
                           - desired_inventory_coverage * expected_order_rate
                           )/7) 
         - 6 * max(0, 1 + (+ desired_inventory_coverage * expected_order_rate
                           - desired_inventory_coverage * expected_order_rate
                           )/7) 
        )/3
        + max(0, 
              1 + (+ 1 * (2 + 2) 
                   - 1 * (2 + 2)
                )/7
            )
      )

    = max(0, 
        (+ 6 * max(0, 1 + (+ (minimum_order_processing_time + safety_stock coverage) * 2 
                           - (minimum_order_processing_time + safety_stock coverage) * 2
                           )/7) 
         - 6 * max(0, 1 + (+ (minimum_order_processing_time + safety_stock coverage) * 2 
                           - (minimum_order_processing_time + safety_stock coverage) * 2
                           )/7) 
        )/3 
        + max(0, 
              1 + (+ 1 * (2 + 2) 
                   - 1 * (2 + 2)
                )/7
            )

      = max(0, 
        (+ 6 * max(0, 1 + (+ (2 + 2) * 2 
                           - (2 + 2) * 2
                           )/7) 
         - 6 * max(0, 1 + (+ (2 + 2) * 2 
                           - (2 + 2) * 2
                           )/7) 
        )/3 
        + max(0, 
              1 + (+ 1 * (2 + 2) 
                   - 1 * (2 + 2)
                )/7
            )
        )

[hyunjimoon]

Vensim provides the following documents:

Backlog= INTEG (
        Order Rate-Order Fulfillment Rate,
            Order Rate * Target Delivery Delay)
    Units: Widgets
    The firm's backlog of unfilled orders

Expected Order Rate= INTEG (
    Change in Exp Orders,
        Customer Order Rate)
Units: Widgets/Week
The demand forecast is formed by adaptive expectations, using 
        exponential smoothing, a common forecasting technique. The 
        initial forecast is equal to the initial customer order rate.

Inventory = INTEG(Production Rate-Shipment Rate,Desired Inventory)
Units: Widgets
The level of finished goods inventory in the plant. Increased by 
        production and decreased by shipments. Initially set to the 
        desired inventory level.

Work in Process Inventory= INTEG (
    Production Start Rate - Production Rate,
        Desired WIP)
Units: Widgets
WIP inventory accumulates the difference between production 
        starts and completions.

Production Rate Stocked= INTEG (
    Production Rate Stocked Change Rate,
        Production Rate)
Units: Widgets/Week

Production Start Rate Stocked= INTEG (
    Production Start Rate Stocked Change Rate,
        Production Start Rate)
Units: Widgets/Week

Backlog= INTEG (
        Order Rate-Order Fulfillment Rate,
            Order Rate * Target Delivery Delay)
    Units: Widgets
    The firm's backlog of unfilled orders

Expected Order Rate= INTEG (
    Change in Exp Orders,
        Customer Order Rate)
Units: Widgets/Week
The demand forecast is formed by adaptive expectations, using 
        exponential smoothing, a common forecasting technique. The 
        initial forecast is equal to the initial customer order rate.

Inventory = INTEG(Production Rate-Shipment Rate,Desired Inventory)
Units: Widgets
The level of finished goods inventory in the plant. Increased by 
        production and decreased by shipments. Initially set to the 
        desired inventory level.

Work in Process Inventory= INTEG (
    Production Start Rate - Production Rate,
        Desired WIP)
Units: Widgets
WIP inventory accumulates the difference between production 
        starts and completions.

Production Rate Stocked= INTEG (
    Production Rate Stocked Change Rate,
        Production Rate)
Units: Widgets/Week

Production Start Rate Stocked= INTEG (
    Production Start Rate Stocked Change Rate,
        Production Start Rate)
Units: Widgets/Week

Backlog= INTEG (
        Order Rate-Order Fulfillment Rate,
            Order Rate * Target Delivery Delay)
    Units: Widgets
    The firm's backlog of unfilled orders

Expected Order Rate= INTEG (
    Change in Exp Orders,
        Customer Order Rate)
Units: Widgets/Week
The demand forecast is formed by adaptive expectations, using 
        exponential smoothing, a common forecasting technique. The 
        initial forecast is equal to the initial customer order rate.

Inventory = INTEG(Production Rate-Shipment Rate,Desired Inventory)
Units: Widgets
The level of finished goods inventory in the plant. Increased by 
        production and decreased by shipments. Initially set to the 
        desired inventory level.

Work in Process Inventory= INTEG (
    Production Start Rate - Production Rate,
        Desired WIP)
Units: Widgets
WIP inventory accumulates the difference between production 
        starts and completions.

Production Rate Stocked= INTEG (
    Production Rate Stocked Change Rate,
        Production Rate)
Units: Widgets/Week

Production Start Rate Stocked= INTEG (
    Production Start Rate Stocked Change Rate,
        Production Start Rate)
Units: Widgets/Week

[Dashadower]

What's done: Codegen no longer blindly substitutes expressions in place of variable names (this led to problems regarding order of operations) function block builder now uses the topological sort class instead of implementing it inside.

TODO: Don't differentiate between stock initialization variables(__init) and other auxiliary variables that's required to evaluate the init values when sorting. Revert stan vensim model interface to use keyword arguments to pass options instead of a dictionary. @hyunjimoon

[hyunjimoon]

@Dashadower

Revert stan vensim model interface to use keyword arguments to pass options instead of a dictionary

My setting_assumption is as follows:

setting_assumption = {
    "est_param_scalar" : ("inventory_adjustment_time", "minimum_order_processing_time"),
    "ass_param_scalar" : ("inventory_coverage", "manufacturing_cycle_time", "safety_stock_coverage", "time_to_average_order_rate", "wip_adjustment_time"),
    "target_simulated_vector" : ( "production_rate_stocked", "production_start_rate_stocked"),
    "driving_data_vector" : ("customer_order_rate", "process_noise_std_norm_data", "production_start_rate_m_noise_trun_norm_data", "production_rate_m_noise_trun_norm_data"),
    "model_name": "mngInven",
    "integration_times": list(range(1, n_t + 1)),
    "initial_time": 0.0
}

For further data analysis, I need the list of classified parameters e.g. est_param_scalarwhich is ["inventory_adjustment_time", "minimum_order_processing_time"]

If we change as your suggestion, could model.est_param_scalar be used simply instead of model.setting_assumption['est_param_scalar']? Without dictionary, I don't how to contain list of somethings.

[hyunjimoon]

Q1. I don't understand StanModelContext which has six input? as follows

class StanModelContext:
    initial_time: float
    integration_times: Iterable[float]
    stan_data: Dict[str, StanDataEntry] = field(default_factory=dict)
    sample_statements: List[SamplingStatement] = field(default_factory=list)
    exposed_parameters: Set[str] = field(default_factory=set)  # stan variables to be passed to the ODE function
    all_stan_variables: Set[str] = field(default_factory=set)  # set of all stan variables

is initialized with only first two input? Could stan_data, sample_statements be initialized later? https://github.com/Data4DM/BayesSD/blob/master/ContinuousCode/5_BayesCalib/stanify/builders/stan/stan_model.py#L169 only uses initial_time, integration_times.

[hyunjimoon]

failure case log

case 1

Tracking work_in_process_inventory_init,

Work in Process Inventory
= Desired WIP
= Manufacturing Cycle Time*Desired Production
= 6.01 * MAX(0,Expected Order Rate+Adjustment from Inventory)
= 6.01 * MAX(0, Customer Order Rate + (Desired Inventory - Inventory)/Inventory Adjustment Time)

Inventory_init = Desired Inventory

= 6.01 * MAX(0, 10.01 + (Desired Inventory Coverage*Expected Order Rate - Desired Inventory Coverage*Expected Order Rate )/7.01)

= 60

which is 
6.01 * MAX(0, 10.01 + (Desired Inventory Coverage*Expected Order Rate - Desired Inventory Coverage*Expected Order Rate )/7.01)

In stan, ((6.01) * (fmax(0, ((10.01)) + ((((0.5) + (2.01)) * ((10.01))) - ((((0.5) + (2.01)) * ((10.01)))) / (7.01)))));

((10.01)) + ((((0.5) + (2.01)) * ((10.01))) - ((((0.5) + (2.01)) * ((10.01)))) / (7.01))))

which is

6.01 *  (10.01 + ((((0.5) + (2.01)) * ((10.01))) - ((((0.5) + (2.01)) * ((10.01)))) / (7.01)))

The following two should be the same:
 ((0.5 + 2.01) * (10.01)) - (((0.5 + 2.01) * 10.01)) / 7.01

 ((0.5 + 2.01) * (10.01)) - (((0.5 + 2.01) * 10.01)) / 7.01

The inventory and its initial value desired_inventory should be crossed out, but not a-a/7 instead of '(a-a)/7` is happening.

case 2

process_noise__init = 1 * 0.1 * 2 - 0.0625 / 10 / 0.0625 / 10 ^ 0.5; 2 - 0.0625 / 10 / 0.0625 / 10 ^ 0.5's equation is ((2-time_step/correlation_time)/(time_step/correlation_time))^0.5 which is not what it is translated.

Process noise's initial value structure is as below:

case3

Different function codes are generated (sorting order I assume) which is why even with the recompile-prevention code if f.read().rstrip() == function_code.rstrip(): it keeps recompile.

case 4

real inventory__init = (170000) * ((0.5) + (2.01));

from draw2data structure, hopefully, there seems be a way to parenthesize in order.

[hyunjimoon]

case 5

data function is translated twice

// Begin ODE declaration
real dataFunc__process_noise_std_norm_data(real time){
    // DataStructure for variable process_noise_std_norm_data
    real slope;
    real intercept;

    if(time <= 1){
        intercept = -0.41079557111850096;
        slope = -0.45109147217853635 - -0.41079557111850096;
        return intercept + slope * (time - -0.41079557111850096);
    }
 ...
}

real dataFunc__process_noise_std_norm_data(real time){
    // DataStructure for variable process_noise_std_norm_data
    real slope;
    real intercept;

    if(time <= 1){
        intercept = -0.41079557111850096;
        slope = -0.45109147217853635 - -0.41079557111850096;
        return intercept + slope * (time - -0.41079557111850096);
    }

    return 197946;
}

real dataFunc__customer_order_rate(real time){
    // DataStructure for variable customer_order_rate
    real slope;
    real intercept;

    if(time <= 1){
        intercept = 146376;
        slope = 147079 - 146376;
        return intercept + slope * (time - 146376);
    }
    else if(time <= 2){
        intercept = 147079;
        slope = 159336 - 147079;
        return intercept + slope * (time - 147079);
    }
    else if(time <= 3){
        intercept = 159336;
        slope = 163669 - 159336;
        return intercept + slope * (time - 159336);
    }
    else if(time <= 4){
        intercept = 163669;
        slope = 170068 - 163669;
        return intercept + slope * (time - 163669);
    }
    else if(time <= 5){
        intercept = 170068;
        slope = 168663 - 170068;
        return intercept + slope * (time - 170068);
    }
    else if(time <= 6){
        intercept = 168663;
        slope = 169890 - 168663;
        return intercept + slope * (time - 168663);
    }
    else if(time <= 7){
        intercept = 169890;
        slope = 170364 - 169890;
        return intercept + slope * (time - 169890);
    }
    else if(time <= 8){
        intercept = 170364;
        slope = 164617 - 170364;
        return intercept + slope * (time - 170364);
    }
    else if(time <= 9){
        intercept = 164617;
        slope = 173655 - 164617;
        return intercept + slope * (time - 164617);
    }
    else if(time <= 10){
        intercept = 173655;
        slope = 171547 - 173655;
        return intercept + slope * (time - 173655);
    }
    else if(time <= 11){
        intercept = 171547;
        slope = 208838 - 171547;
        return intercept + slope * (time - 171547);
    }
    else if(time <= 12){
        intercept = 208838;
        slope = 153221 - 208838;
        return intercept + slope * (time - 208838);
    }
    else if(time <= 13){
        intercept = 153221;
        slope = 150087 - 153221;
        return intercept + slope * (time - 153221);
    }
    else if(time <= 14){
        intercept = 150087;
        slope = 170439 - 150087;
        return intercept + slope * (time - 150087);
    }
    else if(time <= 15){
        intercept = 170439;
        slope = 176456 - 170439;
        return intercept + slope * (time - 170439);
    }
    else if(time <= 16){
        intercept = 176456;
        slope = 182231 - 176456;
        return intercept + slope * (time - 176456);
    }
    else if(time <= 17){
        intercept = 182231;
        slope = 181535 - 182231;
        return intercept + slope * (time - 182231);
    }
    else if(time <= 18){
        intercept = 181535;
        slope = 183682 - 181535;
        return intercept + slope * (time - 181535);
    }
    else if(time <= 19){
        intercept = 183682;
        slope = 183318 - 183682;
        return intercept + slope * (time - 183682);
    }
    else if(time <= 20){
        intercept = 183318;
        slope = 177406 - 183318;
        return intercept + slope * (time - 183318);
    }
    else if(time <= 21){
        intercept = 177406;
        slope = 182737 - 177406;
        return intercept + slope * (time - 177406);
    }
    else if(time <= 22){
        intercept = 182737;
        slope = 187443 - 182737;
        return intercept + slope * (time - 182737);
    }
    else if(time <= 23){
        intercept = 187443;
        slope = 224540 - 187443;
        return intercept + slope * (time - 187443);
    }
    else if(time <= 24){
        intercept = 224540;
        slope = 161349 - 224540;
        return intercept + slope * (time - 224540);
    }
    else if(time <= 25){
        intercept = 161349;
        slope = 162841 - 161349;
        return intercept + slope * (time - 161349);
    }
    else if(time <= 26){
        intercept = 162841;
        slope = 192319 - 162841;
        return intercept + slope * (time - 162841);
    }
    else if(time <= 27){
        intercept = 192319;
        slope = 189569 - 192319;
        return intercept + slope * (time - 192319);
    }
    else if(time <= 28){
        intercept = 189569;
        slope = 194927 - 189569;
        return intercept + slope * (time - 189569);
    }
    else if(time <= 29){
        intercept = 194927;
        slope = 197946 - 194927;
        return intercept + slope * (time - 194927);
    }
    return 197946;
}

vector vensim_ode_func(real time, vector outcome, real minimum_order_processing_time, real inventory_adjustment_time){
    vector[7] dydt;  // Return vector of the ODE function

    // State variables
    real production_start_rate_stocked = outcome[1];
    real production_rate_stocked = outcome[2];
    real backlog = outcome[3];
    real expected_order_rate = outcome[4];
    real work_in_process_inventory = outcome[5];
    real inventory = outcome[6];
    real process_noise = outcome[7];

    real maximum_shipment_rate = inventory / minimum_order_processing_time;
    real target_delivery_delay = 2.01;
    real desired_shipment_rate = backlog / target_delivery_delay;
    real order_fulfillment_ratio = lookupFunc__table_for_order_fulfillment(maximum_shipment_rate / desired_shipment_rate);
    real correlation_time = 10;
    real safety_stock_coverage = 2.01;
    real reference_throughput = 170000;
    real manufacturing_cycle_time = 6.01;
    real production_rate = fmax(0, 1 + process_noise) * work_in_process_inventory / manufacturing_cycle_time;
    real desired_inventory_coverage = minimum_order_processing_time + safety_stock_coverage;
    real desired_inventory = desired_inventory_coverage * expected_order_rate;
    real adjustment_from_inventory = desired_inventory - inventory / inventory_adjustment_time;
    real desired_production = fmax(0, expected_order_rate + adjustment_from_inventory);
    real wip_adjustment_time = 3.01;
    real desired_wip = manufacturing_cycle_time * desired_production;
    real adjustment_for_wip = desired_wip - work_in_process_inventory / wip_adjustment_time;
    real desired_production_start_rate = desired_production + adjustment_for_wip;
    real production_start_rate = fmax(0, desired_production_start_rate);
    real work_in_process_inventory_dydt = production_start_rate - production_rate;
    real shipment_rate = desired_shipment_rate * order_fulfillment_ratio;
    real inventory_dydt = production_rate - shipment_rate;
    real process_noise_scale = 1;
    real order_fulfillment_rate = shipment_rate;
    real order_rate = dataFunc__customer_order_rate(time);
    real dt_over_tc = 0.00625;
    real white_noise = process_noise_scale * dataFunc__process_noise_std_norm_data(time) * 2 - dt_over_tc / dt_over_tc ^ 0.5;
    real process_noise_change_rate = white_noise - process_noise / correlation_time;
    real time_to_average_order_rate = 8;
    real change_in_exp_orders = dataFunc__customer_order_rate(time) - expected_order_rate / time_to_average_order_rate;
    real expected_order_rate_dydt = change_in_exp_orders;
    real backlog_dydt = order_rate - order_fulfillment_rate;
    real time_step = 0.0625;
    real production_rate_stocked_change_rate = production_rate - production_rate_stocked / time_step;
    real production_start_rate_stocked_change_rate = production_start_rate - production_start_rate_stocked / time_step;
    real production_start_rate_stocked_dydt = production_start_rate_stocked_change_rate - production_start_rate;
    real production_rate_stocked_dydt = production_rate_stocked_change_rate - production_rate;
    real process_noise_dydt = process_noise_change_rate;

    dydt[1] = production_start_rate_stocked_dydt;
    dydt[2] = production_rate_stocked_dydt;
    dydt[3] = backlog_dydt;
    dydt[4] = expected_order_rate_dydt;
    dydt[5] = work_in_process_inventory_dydt;
    dydt[6] = inventory_dydt;
    dydt[7] = process_noise_dydt;

    return dydt;
}