Refactoring stanify: Mismatching initialization stock value
Closed this issue · 6 comments
hyunjimoon commented
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_productionadj_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= 1wip_adjustement_time= 3minimum_order_processing_time= 5manufacturing_cycle_time= 6inventory_adjustment_time= 7time_to_average_order_rate= 8safety_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 commented
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 commented
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 commented
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 commented
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 commented
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 commented
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;
}