Solving differential equation with stanify builder?

[hyunjimoon]

Finding equilibrium is needed as we wish to analyze the effect starting from stable (although I wonder what value it might have if real life is far from stability). We usually solve this by hand from chains of equation by assuming every stock's _dydt = 0.

From the file, the following function returns flow (_dydt of each stock) given parameter value minimum_order_processing_time, time_step, inventory_adjustment_time, process_noise_scale. Could we inverse this function and solve for differential equation solution (i.e. value of 7 stocks)?

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

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

    real manufacturing_cycle_time = 8;
    real safety_stock_coverage = 2;
    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 desired_wip = manufacturing_cycle_time * desired_production;
    real wip_adjustment_time = 2;
    real adjustment_for_wip = (desired_wip - work_in_process_inventory) / wip_adjustment_time;
    real desired_production_start_rate = fmax(0, desired_production + adjustment_for_wip);
    real production_start_rate = fmax(0, desired_production_start_rate);
    real desired_minus_shadow_psr = production_start_rate - production_start_rate_stocked;
    real production_start_rate_stocked_change_rate = desired_minus_shadow_psr / time_step;
    real target_delivery_delay = 2;
    real desired_shipment_rate = backlog / target_delivery_delay;
    real production_rate = work_in_process_inventory / manufacturing_cycle_time * fmax(0, 1 + process_noise);
    real desired_minus_shadow_pr = production_rate - production_rate_stocked;
    real production_rate_stocked_change_rate = desired_minus_shadow_pr / time_step;
    real production_rate_stocked_dydt = production_rate + production_rate_stocked_change_rate;
    real time_to_average_order_rate = 8;
    real change_in_exp_orders = (dataFunc__customer_order_rate(time, time_step) - expected_order_rate) / time_to_average_order_rate;
    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 correlation_time_over_time_step = 100;
    real white_noise = 4.89 * correlation_time_over_time_step ^ 0.5 * dataFunc__process_noise_uniform_driving(time, time_step) * process_noise_scale;
    real shipment_rate = desired_shipment_rate * order_fulfillment_ratio;
    real order_fulfillment_rate = shipment_rate;
    real order_rate = dataFunc__customer_order_rate(time, time_step);
    real backlog_dydt = order_rate - order_fulfillment_rate;
    real white_minus_process = white_noise - process_noise;
    real correlation_time = time_step * correlation_time_over_time_step;
    real process_noise_change_rate = white_minus_process / correlation_time;
    real expected_order_rate_dydt = change_in_exp_orders;
    real production_start_rate_stocked_dydt = production_start_rate + production_start_rate_stocked_change_rate;
    real process_noise_dydt = process_noise_change_rate;
    real work_in_process_inventory_dydt = production_start_rate - production_rate;
    real inventory_dydt = production_rate - shipment_rate;

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

    return dydt;
}

How do you usually compute equilibrium @tomfid?

Could stanify builder supply a equation that express stock value at equilibrium as a function of parameter? @Dashadower
To be clear,

• input: length P vector of parameter value (P = 4)
• output: length Q vector of stock value at equilibrium (Q = 7)

[tomfid]

Usually I don't bother computing equilibrium, though I do typically initialize things locally in equilibrium or along a steady-state growth path. This is mainly of interest for experimentation from a known baseline in modest-sized models. Anything calibrated to data is unlikely to be in equilibrium at the start. Also, equilibrium is usually only interesting for a subsystem - in the inventory model, you'd like to see the orders/WIP/inventory in equilibrium, but if you add a cost function, the NPV is a pure accumulation of a nonnegative quantity, so it can't be in equilibrium (except trivially at 0 throughput).

There's no general guarantee of a unique or stable equilibrium, if one even exists. Usually in a complex setting the easiest thing to do is simply pick one or more starting points and run the model for a long time to see what happens, but that may not work. You can also set up a payoff representing the degree of disequilibrium and minimize that, or iteratively linearize and take eigenvectors. Anything else takes some symbolic or structural analysis of the model.

The pred/prey system is a good example - trivial equilibrium with zero populations, and an unstable equilibrium in the middle of the limit cycle orbit. The 0 equilibrium is not really interesting. The unstable point is interesting dynamically, but it's unlikely to exist in a real dataset because noise will always excite the system away from that point. What you'd really like to know is whether some policy can make that point stable, for example.

[hyunjimoon]

Would uniqueness of solution (stable equilibrium) imply non-degeneracy and vice-versa?

[tomfid]

Assuming you mean non-degeneracy in the sense of freedom from pathological posterior behavior, I'd tend to say no.

Example: suppose the model is competitive, with bistable outcomes (A wins, or B wins). The data will be consistent with one or the other, so even though there are 2 equilibria (and a third unstable one in the middle), the likelihood will have 1 maximum.

I'm pretty sure examples of models with 1 equilibrium but multiple maxima also exist, but I can't immediately think of a realistic example.