Data4DM / BayesSD

Data for Decision, Affordable Analytics for All
8 stars 0 forks source link

Check stanify implementation is downsampling #59

Closed hyunjimoon closed 5 months ago

hyunjimoon commented 1 year ago

real production_start_rate_stocked_change_rate = production_start_rate - production_start_rate_stocked / time_step; from function block of inventory model with process noise: is this evaluated 20 vs 20 * (1/0.0625) times?

hyunjimoon commented 1 year ago

As function block is only evaluated for 1,..,20 times, n_t should be (final time - initial time)/ time_step. In this sense, Tc over dt should be the main parameter.

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

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

    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 safety_stock_coverage = 2.01;
    real desired_inventory_coverage = minimum_order_processing_time + safety_stock_coverage;
    real wip_adjustment_time = 3.01;
    real manufacturing_cycle_time = 6.01;
    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 desired_wip = manufacturing_cycle_time * desired_production;
    real adjustment_for_wip = (desired_wip - work_in_process_inventory) / wip_adjustment_time;
    real process_noise_scale = 1;
    real target_delivery_delay = 2.01;
    real desired_shipment_rate = backlog / target_delivery_delay;
    real process_noise_std_norm_data = 0.1;
    real correlation_time = 10;
    real maximum_shipment_rate = inventory / minimum_order_processing_time;
    real order_fulfillment_ratio = lookupFunc__table_for_order_fulfillment(maximum_shipment_rate / desired_shipment_rate);
    real shipment_rate = desired_shipment_rate * order_fulfillment_ratio;
    real production_rate = fmax(0, 1 + process_noise) * work_in_process_inventory / manufacturing_cycle_time;
    real time_step = 0.0625;
    real production_rate_stocked_change_rate = (production_rate - production_rate_stocked) / time_step;
    real desired_production_start_rate = desired_production + adjustment_for_wip;
    real production_start_rate = fmax(0, desired_production_start_rate);
    real production_start_rate_stocked_change_rate = production_start_rate - production_start_rate_stocked / time_step;
    real reference_throughput = 170000;
    real work_in_process_inventory_dydt = production_start_rate - production_rate;
    real order_fulfillment_rate = shipment_rate;
    real production_rate_stocked_dydt = production_rate_stocked_change_rate - production_rate;
    real dt_tc = time_step / correlation_time;
    real white_noise = process_noise_scale * process_noise_std_norm_data * (2 - dt_tc / dt_tc) ^ 0.5;
    real process_noise_change_rate = (white_noise - process_noise) / correlation_time;
    real process_noise_dydt = process_noise_change_rate;
    real inventory_dydt = production_rate - shipment_rate;
    real order_rate = dataFunc__customer_order_rate(time);
    real expected_order_rate_dydt = change_in_exp_orders;
    real backlog_dydt = order_rate - order_fulfillment_rate;
    real production_start_rate_stocked_dydt = production_start_rate_stocked_change_rate - production_start_rate;

    dydt[1] = work_in_process_inventory_dydt;
    dydt[2] = expected_order_rate_dydt;
    dydt[3] = production_rate_stocked_dydt;
    dydt[4] = process_noise_dydt;
    dydt[5] = backlog_dydt;
    dydt[6] = inventory_dydt;
    dydt[7] = production_start_rate_stocked_dydt;
    print(production_rate);
    return dydt;
}