In this blog post, I am going to show you just how easy it is to set up and solve a stochastic programming problem. For the example, I am going to use our Python interface that integrates seamlessly with widely used Python packages such as Numpy and Pandas. Let us begin by importing the needed packages into our work space.
import pandas, quasar
Data and Parameters¶
The decision problem that guides this example is a simple storage optimization model. The decision-maker must decide when and how much to inject and widthdraw from a storage. What makes the decision problem stochastic is that prices follow some random process and any decision taken now effects future opportunities to inject and withdraw as well as how profitable the storage will be. We are planning our storage over next 12 time periods (think of months) and need to account for capacity limits as well as possibly losses.
time_horizon = 12
loss = 0.1
max_fill = 100
max_rate = 25
Model Formulation¶
What is special about QUASAR is that you can formulate any stochastic programming problem just like it was a linear program, but unlike with conventional solvers, any model parameter can be a random variable. All it takes is using the keyword rand along with a string identifier that must be unique for a given stage.
# decision model
model = quasar.DecisionProblem()
# starting state
s = {-1 : 25}
# multiple stages
for t in range(time_horizon):
# decision variables
s[t] = model.add_variable(t, "storage", ub=max_fill)
c = model.add_variable(t, "charge", ub=max_rate)
r = model.add_variable(t, "release", ub=max_rate)
# constraints
s[t] == s[t-1] + (1-loss)*c - r
# objective with random variable
model.obj += quasar.rand("price") * (r - c)
Time Series Data¶
In the above example, prices are labeled as random variables. Now we need to define a stochastic process that describes how prices evolve ofer time. If you happen to have a sample of price scenarios, for example from another model or from real data, simply export these scenarios in Pandas and turn them into time series sample.
df = pandas.read_csv("prices.csv", index_col=[0,1], header=0, names=["price"])
ts = quasar.TimeSeriesSample.from_df(df)
Start the Solver¶
QUASAR's solvers can take it from here. Just pass the model and the time series sample, and the dynamic optimizer will automatically decompose the problem, generate a scenario lattice, and solve it until the solution has converged. This 12-stage problem does not even take a second to be solved.
# optimize problem
opt = quasar.DynamicOptimizer(model, ts)
opt.solve()
opt.summary.tail()
Inspect the Solution¶
The solution of a stochastic program is a decision policy which is itself stochastic (or, more precisely, a function of the scenario as it unfolds over time). The best way to inspect the decision policy is by simulating the optimal policy. In this way, we gain insights into the distribution of decision variables, objective values, and their covariance with the random process. For example, we can use Pandas to receive summary statistics of the simulated objective values at every stage.
# simulate decision process
sim = opt.policy.simulate(sample_size=1000)
# descriptive statistics
sim.df.rewards.groupby(level="stage").describe()
Or we can visualize the quantiles of possible storage fill levels over the entire time horizon by means of a fan chart.
# visulize decision
ax = sim.plot_decision("storage", kind="fanchart")
Of course, this little example just scratches the surface of what QUASAR can do. If you want to see more examples, check out some of our Colab notebooks, which allow you to run QUASAR right in your browser. 