# Stochastic Process Optimization of an Integrated Biorefinery

Dahui Liu1, Sandra Eksioglu1, Mohammad Roni2 1University of Arkansas,
Fayetteville, AR, USA.2Idaho National Laboratory, Idaho Falls, ID,
USA.Manuscript received May 1, 2021.

###### Abstract

Planning of biorefinery operations is complicated by the stochastic nature of
physical and chemical characteristics of biomass feedstock, such as, moisture
level and carbohydrate content. Biomass characteristics affect the performance
of the equipment which feed the reactor and the efficiency of the conversion
process in a biorefinery. We propose a stochastic optimization model to
identify a blend of feedstocks, inventory levels, and operating conditions of
equipment to ensure a continuous flowing of biomass to the reactor while
meeting the requirements of the biochemical conversion process. We propose a
sample average approximation (SAA) of the model, and develop an efficient
algorithm to solve the SAA model. A feedstock preprocessing process consists
of two-stage grinding and pelleting is used to develop a case study. Extensive
numerical analysis are conducted which lead to a number of observations. Our
main observation is that sequencing bales based on moisture level and
carbohydrate content leads to robust solutions that improve processing time
and processing rate of the reactor. We provide a number of managerial insights
that facilitate the implementation of the model proposed.

###### Index Terms:

Production Control, Biomass Processing System, Stochastic Optimization, Sample
Average Approximation, Sequencing, System Reliability

Note for practitioners: This paper is motivated by the challenges faced in the
bioenergy industry. The focus of this paper is on plants which use the
biochemical conversion process to generate liquid fuels. It has been observed
that variations in biomass characteristics, such as moisture content, cause
variations in feeding of the system which lead to under-utilization of
equipment. A requirement of biochemical conversion process is to maintain the
carbohydrate content of biomass processed by the reactor, larger than 59.1%.
We propose a model that identifies the inventory levels and operating
conditions of equipment to ensure a continuous flowing of biomass to the
reactor. The goal is to improve equipment utilization while satisfying the
requirements of the conversion process. The model is tested using real-life
data. We found out that by sequencing bales based on moisture level and
carbohydrate content, a plant can reduce variability in the system leading to
higher processing rates of the reactor.

##  I Introduction

When the Renewable Fuel Standard (RFS) passed in 2007, it set annual target
production amounts for cellulosic biofuel which culminated to 10.5 Billion
Gallons per Year (bgy) by 2020 [2]. These annual targets for cellulosic
biofuel were not met, what led the Environmental Protection Agency (EPA) to
reduce the mandates. A 2016 study by US Department of Energy (DOE) indicates
that US has the necessary resources (biomass) to produce up to 50 bgy of
biofuels by 2040 without jeopardizing the supply of food and animal feed [29].
However, the cost of generating cellulosic biofuels is not competitive to that
of fossil fuels.

Based on [2], the shortfall in meeting the RFS is due to “lack of private
investment, logistical challenges, technology setbacks, and uneven support
from the federal government.” A recent report from US DOE finds “bulk solids
handling and material flows through the system” to be a critical component to
achieve the design throughput of the conversion processes [1]. The flow of
materials is impacted by the variability of biomass characteristics which is
observed when bales of different types of feedstock, with different moisture
level and carbohydrate content, are processed by the same equipment (i.e.,
grinder). These variations negatively affect the integration of biomass
feeding and conversion processes and cause low/unreliable on-stream time of
equipment and low utilization efficiency of the reactor. Indeed, it was
reported by Bell [7] that inefficiencies of handling solid materials are the
reason why biorefineries typically operate at 50% of their capacity in the
first year. Commercial scale cellulosic ethanol producers, such as, Beta
Renewables (the 1st cellulosic ethanol plant in the world), Abengoa,
DowDuPont, etc. have suspended their operations or filed for bankruptcy
because cellulosic biofuel cannot compete with the low price of petroleum
fuel.

The objective of the proposed research is to develop operational strategies
which streamline processes within a biorefinery to ensure a continuous flow of
biomass to the reactor, while meeting the requirements of biochemical
conversion process. This will lead to shorter processing times, higher
utilization of the reactor, and lower processing costs. To accomplish this, we
propose a mathematical model that identifies the infeed rate of biomass to the
system, processing speed of equipment and inventory level to ensure a
continuous feeding of the reactor. It is assumed that different biomass
feedstocks are used in the system, each with different moisture and
carbohydrate contents. Moisture level impacts the infeed rate and processing
speed of equipment since biomass of high moisture level is processed slow to
avoid equipment clogging. Carbohydrate content depends on the feedstock type.
The biochemical conversion process requires that the carbohydrate content of
biomass should be at least 59.1%. Since, both, moisture and carbohydrate
content are random variables, we propose a stochastic optimization model.

Based on our numerical analysis, sequencing bales based on their moisture and
carbohydrate contents reduces variations of the characteristics of biomass
flowing in the system, which could positively impact system’s performance.
Integrating a bale sequencing model within the proposed stochastic
optimization model would make the model challenging to solve. Therefore, we
evaluate different sequences numerically. Via our numerical analysis we
identify characteristics of good sequences and develop a few bale sequencing
rules that are easy to implement. We also provide a number of managerial
insights to facilitate large-scale implementations of the models proposed.

We develop a case study using data collected at the Process Development Unit
(PDU) of Idaho National Laboratory (INL). PDU is a full-size, fully-integrated
feedstock preprocessing system. A number of INL’s industry partners use the
PDU to test different grinding, pelleting, mechanical and chemical separation
options, etc. during design and scale-up of a biorefinery. Figure 1 presents
the flowchart of processes of PDU. Current processes include drag chain
conveyor 1 (DC(1)) to pelleting mill (P). The proposed model also considers a
separation process, which separates biomass based on particle size. Large
particles are processed further. The proposed model identifies a blend of
biomass to meet biochemical conversion specifications as designed by NREL
[13]. These specifications are a total structural carbohydrate of 59.1%, and a
moisture content of 20%.

In summary, these are the main contributions of this research: ($i$) a
stochastic optimization model is developed that leads to robust process
control solutions that meet biochemical process requirements and maintain a
continuous flow of biomass to the reactor; ($ii$) a real-life case study is
developed using historical data from INL’s PDU, and other reliable data
sources; ($iii$) a number of managerial insights are identified that have the
potential to positively impact the performance of biomass processing systems.

The reminder of the paper is organized as follows. In Section II we briefly
review the literature. The proposed stochastic optimization model is
introduced in Section III. In Section IV we present the algorithms developed.
In Section V we present our experimental setup. The results of our experiment
and observations we make are summarized in Section VI. In Section VII we
provide managerial insights about the implementation of the proposed model.

##  II Literature Review

The proposed work aligns with two main streams of literature within the broad
area of supply chain optimization. First, the proposed model is an extension
of the blending problem which identifies minimum cost biomass blends to meet
process requirements of conversion process. Second, the proposed work
contributes to the existing supply chain reliability literature. Below, we
provide details of the literature and highlight our contributions.

###  II-A Biomass Blending

The _blending model_ is a well-known model in operations research [49]. It
identifies a blend of resources to meet demand (or some system requirement) to
minimize costs (or to maximize profits). The model is used in a number of
application in different industries. For example, the model is used in the
agricultural sector to adjust the quality (grade) of grains sold in the grain
market [19]. The model is used in the energy sector to identify coal blends to
meet meet emission criteria [42, 44, 20, 48]. Most recently, a number of
researchers are investigating the opportunities that exist to reduce the cost
of producing biofuels by using biomass blends in biorefineries [41, 43, 15].

The research on biomass blending is motivated by the fact that biomass
feedstocks differ based on their physical characteristics, chemical content,
and selling price. There is an opportunity to reduce the cost of producing
bioenergy by using a blend of feedstocks that meet process requirements at the
minimum cost. However, using biomass blends in a biorefinery is challenging
for the following reasons ($i$) processes should be re-evaluated and re-
adjusted to handle a blend of products with different physical properties,
which impact biomass flowability in the system; ($ii$) the efficiency of the
conversion process is affected by the chemical composition of the biomass
blend; ($iii$) business relationships should be built and maintained with a
number of suppliers to obtain feedstocks; etc. In response to these
challenges, three main streams of research have been developed that focus on
evaluating conversion process [12, 39, 50], evaluating the economic
feasibility and risks within the biorefinery [46, 36, 28, 26], and evaluating
the supply chain costs and risks of using biomass blends [22, 41, 33, 40, 10,
27].

The proposed research contributes to the literature focused on evaluating the
economic feasibility of biomass blend within a biorefinery. Different from the
literature, the proposed work demonstrates the impact of biomass blending on
flow, reactor’s utilization, inventory level, etc. We propose strategies to
sequence bales which ensure a continuous flow of biomass to the reactor.

###  II-B Systems Reliability

A recent work by [11] provides a review of the reliability literature and
classifies existing models into ($i$) redundancy allocation, ($ii$)
reliability allocation, ($iii$) reliability-redundancy allocation, and ($iv$)
assignment and sequencing. Our proposed model contributes to redundancy
allocation and sequencing streams of research since it identifies the amount
of inventory to ensure continuous feeding of the reactor. In the numerical
analysis, we investigate the impact of bale sequencing in the performance of
the system.

The model proposed falls under ($i$) redundancy allocation and ($iv$)
assignment and sequencing streams of research for the following reasons. Our
model identifies buffer size and location, that add redundancies to the
system. This inventory contributes to maintaining a continuous flow of biomass
to the reactor. We have evaluated numerically the impact of bale sequencing
based on moisture level and carbohydrate content to the performance of the
system. Other researchers use extensions of the redundancy allocation model to
design processes in a bioenergy plant [6, 4, 5, 18, 9]. Different from these
works, our model considers the impact of bale sequencing, and biomass blending
to processes in a bioenergy plant.

The model we propose is an extension of the work by [25] that identifies a
process control strategy to streamline preprocessing of switchgrass in a
biorefinery. Different from [25], we propose a stochastic optimization model
that considers the impacts of stochastic carbohydrate and moisture content to
the performance of the system. Additionally, our proposed model considers a
blend of biomass feedstocks.

Finally, the model we propose is an extension of the flow shop scheduling
problem that identifies a sequence of (bales) jobs processed that minimizes
makespan [16, 17, 51, 32, 35]. The order of operations is the same for all
(bales) jobs, but the processing time of each (bale) job in each operation
differs based on the characteristics of the (bale) job. In our model, we group
bales based on moisture level (low, medium, high) and carbohydrate content
($<59.1\%$ or $\geq 59.1\%$). Since this grouping leads to only a few
different types of bales, we evaluate sequences numerically.

##  III Model Formulation

###  III-A Problem Description

Figure 1 presents a flowchart of processes at the PDU. The current processes
include DC(1), which feeds biomass bales to the system, to P where biomass is
pelleted. DC(1) moves biomass to Grinder 1 (G(1)). A separation process after
G(1) separates biomass based on particle size. Large particles are moved via
DC(2) and DC(3) to G(2) for further processing, and from G(2) processed
biomass is moved to metering bin (M). Small particles are moved via SC(1) to
M. Via SC(5), biomass is moved from M to pelleting mill (P).

The proposed flowchart includes additional storage (S(1) to S($n$)) and
conveyors (DC(4) to DC(2$n$ +3)). Each additional storage unit inventories
pellets of a particular biomass feedstock. Each feedstock has specific
carbohydrate and ash contents. A mix of pellets, which meet biochemical
conversion process requirements, is fed to the reactor (R) in every time
period.

Figure 1: Proposed Flowchart of Processes at PDU.

###  III-B A Mixed-Integer Nonlinear Model

Let us assume that $n_{sb}$ bales of feedstock $b\in\mathcal{B}$, with
moisture level $s\in\mathcal{S}$, weighting $m_{sb}$ dry tons (dt) each, are
to be processed at the PDU. Let $t(\in\mathcal{T})$ represent a time period
within the planning horizon. The objective is to identify how much biomass to
process in every time period and how much biomass to carry in the inventory in
order to minimize the time it takes to process every bale. Let $\tilde{Z}_{t}$
be a binary variable that takes the value 1 if the reactor is up and running
in period $t$, and takes the value 0 otherwise. The objective of our model is

$\displaystyle\min$ $\displaystyle\sum_{t\in\mathcal{T}}\tilde{Z}_{t}.$ (1)

Biomass processing is impacted by the capacity of equipment used. Let
$u_{isb}$ represent the capacity of equipment $i(\in\mathcal{I})$ in terms of
dt of biomass processed per time period. The capacity depends on the type of
biomass processed and its moisture level. For example, the PDU operator lowers
the rotational speed of grinders when processing highly moist biomass in order
to avoid equipment clogging. Thus, the processing capacity of an equipment is
higher for dry biomass then for moist biomass. Additionally, the capacity of
storage bins ($\mathcal{I^{M}}\cup\mathcal{I^{B}}$) is expressed in terms of
mass ($m_{i}$) and volume ($v_{i}$), because the density of processed biomass
($d_{isb}$) impacts the space it takes in storage. Biomass feedstocks differ
based on density, which impacts the amount of biomass that can be stored in
the same storage space, thus, equipment capacity.

Let $X_{isbt}$ represent the amount of biomass flowing from equipment $i$, and
let $M_{isbt}$ represent the amount of biomass in inventory. (2) to (5)
represent the capacity constraints of processing, transportation and storage
equipment.

$\displaystyle X_{1sbt}\leq{u}_{1sb}{Z}_{sbt},$ $\displaystyle\forall
s\in\mathcal{S},b\in\mathcal{B},t\in\mathcal{T},$ (2) $\displaystyle
X_{isbt}\leq{u}_{isb},$ $\displaystyle\forall
i\in\mathcal{I},s\in\mathcal{S},b\in\mathcal{B},t\in\mathcal{T},$ (3)
$\displaystyle\sum_{s\in\mathcal{S}}\sum_{b\in\mathcal{B}}M_{isbt}\leq{m}_{i},$
$\displaystyle\forall i\in\mathcal{I^{M}}\cup\mathcal{I^{B}},t\in\mathcal{T},$
(4)
$\displaystyle\sum_{s\in\mathcal{S}}\sum_{b\in\mathcal{B}}\frac{M_{isbt}}{{d}_{isb}}\leq{v}_{i},$
$\displaystyle\forall i\in\mathcal{I^{M}}\cup\mathcal{I^{B}},t\in\mathcal{T}.$
(5)

In constraint (2), ${Z}_{sbt}$ is binary decision variable which takes the
value 1 if biomass feedstock $b$ is processed in period $t$, and takes the
value 0 otherwise. The value of these variables is determined by constraint
(6). When ${Z}_{sbt}=0,$ the feeding of the system ends, what eventually ends
reactor’s processing.

$\displaystyle\sum_{s^{\prime}\in\mathcal{S}}\sum_{b^{\prime}\in\mathcal{B}}\sum_{t^{\prime}=t}^{t+p_{sb}-1}{Z}_{s^{\prime}b^{\prime}t^{\prime}}-\sum_{t^{\prime}=t}^{t+p_{sb}-1}{Z}_{sbt^{\prime}}$
$\displaystyle\leq M(1-y_{sbt}),$ $\displaystyle\forall
s\in\mathcal{S},b\in\mathcal{B},t\in\mathcal{T}.$ (6)

In constraint (6), $y_{sbt}$ is a parameter which takes the value 1 if a bale
with characteristics $b$ and $s$ is processed in period $t$, and takes the
value zero otherwise. The values of this parameter represent the sequence of
bales fed to the system. In this model we consider this to be a problem input.
Nevertheless, mathematical models can be used to sequence the bales.

The following constraints mimic how the system operates. Constraint (7)
determines the amount of biomass fed to the system in period $t$. This amount
depends on the infeed rate of the system (in meters per time unit), $V_{sbt}$,
and the amount of biomass in (each meter of) a bale, $c_{sb}$. To estimate
$c_{sb}(=w\times h\times d_{sb})$ we use the dimensions of a bale ($w,h$) and
its density ($d_{sb}$). Constraint (8) ensures that the amount of biomass fed
to the system during the planning horizon does not surpass the amount
available. Constraints (9) and (10) present the relationship between the
infeed rate and the time it takes to process a bale of length $l$. Constraint
(11) determines the values of $\tilde{Z}_{t},$ and helps us determine the
minimum processing time of the reactor.

$X_{1sbt}\leq c_{sb}V_{sbt},\forall
s\in\mathcal{S},b\in\mathcal{B},t\in\mathcal{T},$ (7)
$\sum_{t\in\mathcal{T}}X_{1sbt}\leq m_{sb}n_{sb},\forall
s\in\mathcal{S},b\in\mathcal{B},$ (8)

$\sum_{t^{\prime}=t}^{t+p_{sb}-1}V_{sbt^{\prime}}+{M}(1-y_{sbt})\geq l,\forall
s\in\mathcal{S},b\in\mathcal{B},t\in\mathcal{T},$ (9)
$\sum_{t^{\prime}=t}^{t+p_{sb}-1}V_{sbt^{\prime}}-{M}(1-y_{sbt})\leq l,\forall
s\in\mathcal{S},b\in\mathcal{B},t\in\mathcal{T},$ (10)

$\tilde{Z}_{t}\leq\tilde{Z}_{t-1}\forall t\in\mathcal{T}\backslash\\{1\\}.$
(11)

Constraints (12) to (15) represent the flow balance constraints. In (12), the
amount of flow from grinder $i$ is reduced based on dry matter loss, $\mu_{i}$
(%). Constraints (13) and (14) model the separation process. Here,
$\theta_{sb}$ is the bypass ratio of this system.

$X_{isbt}=\sum_{p\in\mathcal{I}_{i}}(1-\mu_{i})X_{psbt},\forall
i\in\mathcal{I^{G}},s\in\mathcal{S},b\in\mathcal{B},t\in\mathcal{T},$ (12)

$\displaystyle
X_{ksbt}=\sum_{p\in\mathcal{I}_{i}}\left(1-\theta_{sb}\right)X_{psbt},k=\text{DC(2)},\forall
s\in\mathcal{S},b\in\mathcal{B},$ $\displaystyle t\in\mathcal{T},$ (13)

$X_{ksbt}=\sum_{p\in\mathcal{I}_{i}}\theta_{sb}X_{psbt},k=\text{SC(1)},\forall
s\in\mathcal{S},b\in\mathcal{B},t\in\mathcal{T},$ (14)

$\displaystyle X_{isbt}=\sum_{p\in\mathcal{I}_{i}}X_{psbt}\forall
i\in\mathcal{I}\backslash\\{\text{DC(2),
SC(1)},\mathcal{I^{M}},\mathcal{I^{B}},\mathcal{I^{G}}\\},$ $\displaystyle
s\in\mathcal{S},b\in\mathcal{B},t\in\mathcal{T}.$ (15)

Constraint (16) presents the flow conservation constraint in storage bins.
Constraints (17) and (18) show that in initial and ending inventory are zero.

$\displaystyle
M_{isbt}=M_{isbt-1}+\sum_{p\in\mathcal{I}_{i}}X_{psbt}-X_{isbt}$
$\displaystyle\forall i\in\mathcal{I^{M}\cup I^{B}},$ $\displaystyle
s\in\mathcal{S},b\in\mathcal{B},t\in\mathcal{T},$ (16) $\displaystyle
M_{isb0}=0,\forall i\in\mathcal{I^{M}\cup
I^{B}},s\in\mathcal{S},b\in\mathcal{B},$ (17) $\displaystyle
M_{isbT}=0,\forall i\in\mathcal{I^{M}\cup
I^{B}},s\in\mathcal{S},b\in\mathcal{B}.$ (18)

Constraints (19) to (21) model reactor’s reliability. These constraints ensure
that the feeding of the reactor is within a certain bound of its maximum
capacity. Let $U$ represent the maximum feeding rate of the reactor.
Constraints (19) and (20) ensure that the feeding of the reactor is within
$[\underline{q}U,U]$ (for $0\leq\underline{q}<1$) in every time period.
Additionally, (21) ensures that, on the average, the feeding of the reactor is
at least $qU$ for $q\in(\underline{q},1]$. We multiply the right-hand-side of
(20) and (21) by $\tilde{Z}_{t}$ to estimate the amount of biomass processed
during the time that the reactor is up and running (i.e., $\tilde{Z}_{t}=1$).

$\displaystyle\sum_{i\in\mathcal{I}_{r}}\sum_{s\in\mathcal{S}}\sum_{b\in\mathcal{B}}X_{isbt},$
$\displaystyle\leq{U},$ $\displaystyle\forall t\in\mathcal{T},$ (19)
$\displaystyle\sum_{i\in\mathcal{I}_{r}}\sum_{s\in\mathcal{S}}\sum_{b\in\mathcal{B}}X_{isbt}$
$\displaystyle\geq\underline{q}{U}\tilde{Z}_{t},$ $\displaystyle\forall
t\in\mathcal{T},$ (20)
$\displaystyle\sum_{i\in\mathcal{I}_{r}}\sum_{s\in\mathcal{S}}\sum_{b\in\mathcal{B}}\sum_{t\in\mathcal{T}}X_{isbt}$
$\displaystyle\geq{q}U\sum_{t\in\mathcal{T}}\tilde{Z}_{t}.$ (21)

Let $\mathcal{X}$ represent the feasible region of the optimization model
described via constraints (2) to (21), and ${\bf x}\in\mathcal{X}$ a solution
in this region. Let $f({\bf x})$ represent the objective (1) as a function of
${\bf x}$. The following is a succinct formulation of the optimization model
described in this Section

$\min:f({\bf x})\mbox{ s.t. }{\bf x}\in\mathcal{X}.$

Let us refer to this as model formulation ($P$).

###  III-C A Linear Reformulation of Model ($P$)

In ($P$), the maximum feeding rate of the reactor, $U$, is a decision
variable. We do not set this variable to a predetermined value in order to
allow the model identify the maximum feeding rate of the reactor for given the
available storage space, the quality of biomass and the processing speed of
equipment.

In constraints (20) and (21) of ($P$) the terms ${U}\tilde{Z}_{t}$ are
bilinear, thus, ($P$) is a mixed-integer nonlinear program (MINLP). We use the
McCormick relaxation [31] of these bilinear terms. The corresponding
relaxation is a mixed integer linear program (MILP). Furthermore, since
$\tilde{Z}_{t}$ are binary variables, this relaxation is tight [3].

Let ($\bar{P}$) be the new MILP. In ($\bar{P}$), $W_{t}={U}\tilde{Z}_{t}$.
Thus, the right-hand-side of (20) and (21) are updated to $\underline{q}W_{t}$
and $q\sum_{t\in\mathcal{T}}W_{t}$, respectively. Additionally, the following
constraints are added to the constraint set.

$\displaystyle W_{t}$ $\displaystyle\geq\underline{U}\tilde{Z}_{t},\hskip
79.49744pt\forall t\in\mathcal{T},$ (22) $\displaystyle W_{t}$
$\displaystyle\leq\overline{U}\tilde{Z}_{t},\hskip 79.49744pt\forall
t\in\mathcal{T},$ (23) $\displaystyle W_{t}$
$\displaystyle\geq\overline{U}\tilde{Z}_{t}+U-\overline{U},\hskip
36.135pt\forall t\in\mathcal{T},$ (24) $\displaystyle W_{t}$
$\displaystyle\geq\underline{U}\tilde{Z}_{t}+U-\underline{U},\hskip
36.135pt\forall t\in\mathcal{T}.$ (25)

In constraints (22) to (25), $\underline{U},\overline{U}$ are fixed and known
lower and upper bounds of $U.$ Let $\bar{\mathcal{X}}$ represent the feasible
region of ($\bar{P}$).

###  III-D A Chance Constraints Model

Blending of biomass feedstock with different physical and chemical
characteristics provides an opportunity to meet the requirements of the
conversion platform at minimum cost. For example, a report by the National
Renewable Energy Laboratory (NREL) [13] sets the design specifications for the
biochemical conversion. These specification recommend that the total
structural carbohydrate of the feedstock used should be at least 59.1%.

Figure 2: Distribution of Carbohydrates (% glucan plus % xylan) for Biomass
Feedstocks [37].

A total of 2,000 samples of biomass was collected and analyzed during 2008 and
2014 by [37]. A summary of the carbohydrate contents of these samples is
presented in Figure 2. Based on this data, switchgrass and mischanthus have
higher carbohydrate contents as compared to other feedstocks analyzed. The
average carbohydrate contents for these feedstocks are close to the 59.1%
threshold, which means that about 50% of the samples collected do not meet the
threshold. Although biorefineries would prefer to use biomass feedstocks with
carbohydrate contents higher then this threshold all the time, achieving this
may be impractical and costly. Thus, in our model setting we assume that a
biorefinery strives to meet the carbohydrate content threshold most (i.e., 80%
or 90%) of the time. To model this, we use chance constraints.

Let $f_{b}$ be a stochastic parameter which represents the carbohydrate
content of biomass feedstock $b$. Let $f^{*}$ represent the threshold level
that the biorefinery strives to achieve, and let $\gamma$ represent the risk
level associated with achieving this threshold. Constraint (26) is the chance
constraint which indicates that, during the planning horizon $\mathcal{T}$,
the carbohydrate threshold requirements are met at least $(1-\gamma)\%$ of the
time.

$\displaystyle
Pr\left(\sum_{i\in\mathcal{I}_{r}}\sum_{s\in\mathcal{S}}\sum_{b\in\mathcal{B}}\sum_{t\in\mathcal{T}}(f_{b}-f^{*})X_{isbt}\geq
0\right)$ $\displaystyle\geq 1-\gamma.$ (26)

Let $H({\bf
x},f)=\sum_{i\in\mathcal{I}_{r}}\sum_{s\in\mathcal{S}}\sum_{b\in\mathcal{B}}\sum_{t\in\mathcal{T}}(f^{*}-f_{b})X_{isbt}.$
The following is a succinct formulation of the chance constraint optimization
model, which we refer to as model ($Q$).

$\min_{\bf x\in\mathcal{X}}:f({\bf x})\mbox{ s.t. }Pr\left(H({\bf x},f)\geq
0\right)\leq\gamma.$

###  III-E A Sample Average Approximation Model

Chance constraint models, such as ($Q$), are difficult to solve for two main
reasons [8]. First, calculating the probability of meeting a certain
constraint is often challenging due to multidimensional integration. Second,
the feasibility region defined by the chance constraints may not be convex
[24]. There are two main approaches solve chance constraints models. Both
approaches approximate the probabilistic constraints. One approach discretizes
the probability distribution function of the random parameters, and solves the
corresponding combinatorial problem [14, 30]. The other approach develops
convex approximations of the probabilistic constraint [34]. We follow the
first approach to solve ($Q$).

In model ($Q$), the probability density function of $f_{b}$ does not follow a
known distribution (see Figure 2), thus, calculating $Pr\left(H({\bf x},f)\geq
0\right)$ is challenging. We use the Sample Average Approximation (SAA) method
to approximate the probability density function of $f_{b}$ via an empirical
density function. Let $\hat{H}_{N}({\bf
x},f)=N^{-1}\sum_{n=1}^{N}\Delta(H_{n}({\bf x},f_{n}))$ denote an empirical
measure of the probability distribution of $H({\bf x},f).$ Here,
$f_{n}=\\{f_{bn}|b\in\mathcal{B}\\}$ (for $n=1,\ldots,N$) are $N$ iid samples
of $f_{b}$; $\Delta()$ is a measure of the probability mass function value;
and $N^{-1}$ is the probability assigned to observations of $f_{b}$.

We can approximate $H({\bf x},f)$ using the empirical measure
$\hat{H}_{N}({\bf x},f)$ as follows.

$\displaystyle\bar{H}_{N}({\bf
x})=\mathbb{E}\bigg{[}\mathbb{1}_{(0,\infty)}\big{(}\hat{H}_{n}({\bf
x},f_{n})\big{)}\bigg{]}=$
$\displaystyle\frac{1}{N}\sum_{n=1}^{N}\mathbb{1}_{(0,\infty)}(\hat{H}_{n}({\bf
x},f_{n})).$ (27)

Where, $\mathbb{1}_{(0,\infty)}(t):=\bigg{\\{}\begin{array}[]{c}1\mbox{ if
}t>0\\\ 0\mbox{ if }t\leq 0\end{array}.$ $\bar{H}_{N}({\bf x})$ returns what
proportion of the time $\hat{H}({\bf x},f)\geq 0.$ The resulting SAA
approximation of ($Q$), referred to as model ($\bar{Q}$), is presented below.

$\min_{\bf x\in\mathcal{X}}:f({\bf x})\mbox{ s.t. }\bar{H}_{N}({\bf
x})\leq\hat{\gamma}.$

Note that $\hat{\gamma}$ is the risk level of ($\bar{Q}$). Based on [30], if
$\hat{\gamma}<\gamma$, then, the probability that a feasible solution of model
($\bar{Q}$) is feasible to the true model ($Q$) approaches 1 as
$N\rightarrow\infty.$

##  IV Solution Approach

###  IV-A An Algorithm to Solve Model ($\bar{Q}$)

The constraint set of ($\bar{Q}$) includes an indicator function. Since
commercial software cannot handle such a function, we reformulate these
constraints by introducing the following continuous variables, $\beta^{-}$ and
$\beta^{+}.$ These variables quantify the violation of constraint
$\bar{H}_{N}({\bf x})\leq\hat{\gamma}$. In the following model we have
approximated this constraint using linear functions.

$\displaystyle\nu^{*}(\alpha)=\min_{\bf x\in\mathcal{X}}:f({\bf
x})+\alpha\sum_{n=1}^{N}\beta_{n}^{-}$ (28) Subject to: $\displaystyle{H}({\bf
x},f_{n})+\beta_{n}^{+}-\beta_{n}^{-}$ $\displaystyle=0$ $\displaystyle\forall
n=1,\ldots,N,$ (29) $\displaystyle\beta_{n}^{-},\beta_{n}^{+}$
$\displaystyle\geq 0$ $\displaystyle\forall n=1,\ldots,N.$ (30)

We refer to this formulation as model ($\hat{Q}_{\alpha}$). Its objective
minimizes the cost of violating the chance constraint. The value of the
parameter $\alpha$ in the objective function determines the penalty for
violating the chance constraint. When the value of $\alpha$ is very large, the
objective sets the value of $\beta^{-}_{n}=0$ for all $n=1,\ldots,N.$
Consequently, $H_{n}({\bf x},f_{n})\leq 0$ for all $n=1,\ldots,N.$ When the
value of $\alpha$ is very small, the objective sets the value of
$\beta^{-}_{n}\geq 0$ which leads to solutions which violate the chance
constraint. Thus, the goal is to find
$0<\alpha\in\big{[}\alpha^{l},\alpha^{u}\big{]}$ which minimizes
$\nu^{*}(\alpha)$ while ensuring that the chance constraints are satisfied in
at least $\lceil(1-\hat{\gamma})N\rceil$ of the samples generated.

$\displaystyle\min_{\alpha}:\nu^{*}(\alpha)$ (31) Subject to: $\displaystyle
0<\alpha\in\big{[}\alpha^{l},\alpha^{u}\big{]},$ (32)
$\displaystyle\sum_{n=1}^{N}\mathbb{1}_{(0,\infty)}\left(\beta_{n}^{-}\right)\geq\lceil(1-\hat{\gamma})N\rceil.$
(33)

We refer to this formulation as model ($\hat{Q}$). The SAA Algorithm described
below solves ($\hat{Q}$) using a binary search approach.

Algorithm 1 SAA Algorithm

Initialization: $0<\alpha^{l}\leq\alpha\leq\alpha^{u}$; $\epsilon=10^{-4},$
$\phi=10^{-4}$; Found = False.

1:while (Found = False) do

2: $\alpha\leftarrow\frac{\alpha^{l}+\alpha^{u}}{2}$

3: Solve ($\hat{Q}$); Let $\hat{\beta}_{n}^{-}$ be the incumbent solution of
$\beta_{n}^{-}$

4: $C\leftarrow 0$

5: for $n=1,\ldots,N$ do

6: if $\beta_{n}^{-}>0$ then

7: $C\leftarrow C+1$

8: end if

9: end for

10: if $C\geq\hat{\gamma}N+\epsilon$ then

11: $\alpha^{l}\leftarrow\alpha$

12: else if $C<\hat{\gamma}N-\epsilon$ then

13: $\alpha^{u}\leftarrow\alpha$

14: end if

15: if $|\alpha-\frac{\alpha^{l}+\alpha^{u}}{2}|\leq\phi$ then

16: Found = True

17: end if

18:end while

19:Return $\alpha$ and the solution of ($\hat{Q}$)

###  IV-B An Algorithm to Solve Model ($Q$)

Let us briefly describe the algorithms we use to generate probabilistic lower
and upper bounds for the true chance constraint model ($Q$). We follow closely
the approach proposed by [30].

Let $\nu^{*}$ and $\hat{\nu}_{N}$ represent the optimal objective function
value of ($Q$) and ($\hat{Q}$), respectively. Since ($\hat{Q}$) is a
relaxation of ($Q$), and for $\hat{\gamma}>\gamma$, $\hat{\nu}_{N}$ is a lower
bound for $\nu^{*}$. Indeed, the probability that $\hat{\nu}_{N}$ is a lower
bound to $\nu^{*}$ approaches 1 as $N\rightarrow\infty.$ If
$\hat{\gamma}<\gamma$, then, the probability that $\hat{\nu}_{N}$ is an upper
bound for $\nu^{*}$ approaches 1 as $N\rightarrow\infty.$ Work by [30]
determines the sample size $N$, the number of replications $M$, and the
following approaches to find high confidence upper/lower bounds for ($Q$).

_Lower Bound Procedure_ : We follow these steps: (1) Select
$\hat{\gamma}>\gamma$, and
$N\geq\frac{1}{2(\hat{\gamma}-\gamma)^{2}}\log{\big{(}\frac{1}{\delta}\big{)}}$.
(2) Solve $M=10$ replications of ($\hat{Q}$). In each replication we use a
different sample of size $N.$ Let $\hat{\nu}_{Nj}$ represent the corresponding
objective function values for $j=1,\ldots,M.$ (3) Find the
$\underline{\nu}=min_{j=1,\ldots,M}\hat{\nu}_{Nj}.$ Then, $\underline{\nu}$
represents the lower bound of ($Q$) with confidence $(1-\delta)$.

_Upper Bound Procedure_ : We follow these steps: (1) Select $N$, $\delta$,
$\hat{\gamma}<\gamma$, $M=10$. (2) Solve ($\hat{Q}$) for each of the $M$
replications, and let ${\bf x}_{jN}$ represent the solution found
($j=1,\ldots,M$). (3) Conduct an _a posteriori_ check to see if the chance
constraint is satisfied, that is, $Pr\left(H({\bf x}_{jN},f)\geq
0\right)\leq\gamma.$ In order to estimate this probability, we sample
$N^{\prime}$ ($>>N$) iid values of the random parameter $f_{b}.$ We use
equation (III-E) to calculate $\bar{H}_{N^{\prime}}({\bf x}_{jN},f)$, which is
an estimate of the $Pr\left(H({\bf x}_{jN},f)\geq 0\right)$. (4) Calculate
$\mathcal{U}_{N^{\prime}}({\bf x}_{jN},f)=\bar{H}_{N^{\prime}}({\bf
x}_{jN},f)+\Phi^{-1}(1-\delta)\sqrt{\frac{\big{(}\bar{H}_{N^{\prime}}({\bf
x}_{jN},f)\big{)}\big{(}1-\bar{H}_{N^{\prime}}({\bf
x}_{jN},f)\big{)}}{N^{\prime}}},$ where $\Phi^{-1}()$ represents the inverse
c.d.f. of $N(0,1)$. If $\mathcal{U}_{N^{\prime}}({\bf x}_{jN},f)\leq\gamma,$
then, ${\bf x}_{jN}$ is feasible to ($Q$) with probability ($1-\delta$). (5)
If every ${\bf x}_{jN}$ ($j=1,\ldots,M$) is feasible, then, STOP; otherwise,
increase $N$ and repeat steps (2) to (5).

The pseudo codes of the lower and upper bound procedures are presented in
Appendix A.

###  IV-C Distributed and Parallel Implementation of the SAA Algorithm

In our computational analysis, we solve $M$ replications of ($\hat{Q}$) to
generate bounds for the true problem ($Q$). Solving ($\hat{Q}$) via the SAA
Algorithm is time consuming because the MILP formulation ($\hat{Q}_{\alpha}$)
is solved
$\big{\lceil}\frac{\ln_{2}(\alpha_{u}-\alpha_{l})}{\ln_{2}(\phi)}\big{\rceil}$
times; and each time with a different value of $\alpha$. Thus, we solve $M$
problems using distributed and parallel computing capabilities of the Plametto
cluster, a high-performance computing resource of Clemson University.

Via distributed computing, the computational problem is separated into several
individual tasks, where each of task is completed on one or more nodes. Unlike
the distributed computing, in parallel computing tasks are further divided
into smaller tasks which are then solved using multi-processors on a single
computer [45]. The literature presents several different parallel computing
approaches, including bit-level, instruction-level, data, and task parallelism
[38]. We use task parallelism, which execute different tasks across multiple
threads at the same time.

In the distributed computing environment of the Palemtto cluster, let $N^{c}$
represent the computing nodes (computers) that communicate with one another to
solve the $M$ different problems (replications), each assigned to a node.
These problems are solved in parallel, that is, each problem is solved
simultaneously and without any intersection or disturbance. Let $P^{t}$
represent the total number of processors per node. In Appendix A we present
the pseudo code of the distributed and parallel computing procedure we
develop.

##  V Case Study and Experimental Setup

###  V-A Case Study

We develop a case study to validate the models developed and derive managerial
insights. The case study considers that a blend of corn stover 2-pass (2-P),
corn stover 3-pass (3-P), switchgrass and miscanthus is processed in a
facility which uses biochemical conversion to produce renewable fuel. Figure 1
presents the flowchart of processes in this facility. This process is expected
to deliver a consistent feedstock which meets biochemical conversion
specifications as designed by the National Renewable Energy Laboratory with a
total structural carbohydrate of 59.1% and a moisture content of 20% [13]. We
utilize historical data collected at INL’s PDU during the recent few years.
Table XVI in Appendix summarizes technical, process related data for
switchgrass.

Figures 3 and 4 represent the empirical distribution of carbohydrate content
for each feedstock. In order to develop these empirical distributions we used
historical collected at INL, data collected at DuPont [23], and data reported
by [37].

Figure 3: Empirical distribution of corn stover 3-P and corn stover 2-P.

Figure 4: Empirical distribution of swichgrass and miscanthus. TABLE I: Input
Data: Biomass Blend

| Nr. of Bales |
---|---|---
Feedstock | Moisture Level | Total | Avg. Carb.
Type | L | M | H |  | Cont. (%)
Corn Stover (3-P) | 6 | 10 | 4 | 20 | 57.4
Corn stover (2-P) | 12 | 20 | 8 | 40 | 60.3
Switchgrass | 3 | 5 | 2 | 10 | 66.6
Miscanthus | 3 | 5 | 2 | 10 | 81.7

Table I summarizes the data we use in the case study. We consider that 80
bales of biomass are to be processed in this system. For each feedstock type
we report the total number of bales, the number of bales per moisture level,
and the average carbohydrate content per bale. Table II presents the sequences
of bales we use. The performance of the system is impacted by the sequence,
thus, we solve 6 problems, each using a different sequence. In these sequences
we use letters H, L, M to represent high, low and medium moisture content
correspondingly. We use letters C2, C3, S and M to represent corn stover 2-P,
stover 3-P, switchgrass and miscanthus correspondingly. The numbers represent
the number of bales. We also use a random sequence, and a short sequence. The
short sequence is presented in the Appendix.

TABLE II: Problem Definitions

Problem | Sequence Pattern
---|---
Nr. | Feedstock-Based | Moisture-Based
1 | 10S-40C2-10M-20C3 | 3L-3M-2H
2 | 10S-20C2-10M-20C2-20C3 | 3L-3M-2H
3 | Short sequence | 3L-3M-2H
4 | Short sequence | 5M-3L-2H
5 | Short sequence | 2H-3L-5M
6 | Random sequence | 3L-3M-2H

###  V-B Experimental Setup

The algorithms proposed are coded in Python [47]. Gurobi Optimizer is used to
solve the MILPs [21]. The parallel implementation of SAA Algorithm was
conducted in the Plametto cluster of Clemson University.

Tables III and IV summarize the results from our computational effort to
identify a sample size $N$ and $\hat{\gamma}$ for which the solution of model
($\hat{Q}$) is feasible for model ($Q$), with high confidence level. We follow
the procedure described in Section IV-B. In these experiments, $M=10,$ and
$\gamma\geq\hat{\gamma}.$ We only report the results of solving problem 3 for
different values of $N$ and $\hat{\gamma}$. The confidence level is
$(1-\delta)*100=99\%.$

Based on our numerical results, when $\gamma>\hat{\gamma}$ (Table III) and $N$
= 400 or $N$ = 600, each of $M$ replications gave a feasible solution. The
corresponding maximum solution risk is smaller then $\gamma.$ When
$\gamma=\hat{\gamma}$ (Table IV) and $N$ = 2,500, each of $M$ replication gave
a feasible solution. The corresponding maximum solution risk is smaller then
$\gamma.$ Notice that the computational efforts to find feasible solutions
when $\gamma=\hat{\gamma}$ are higher then in the case when
$\gamma>\hat{\gamma}$. The quality of the optimal solutions for
$\gamma=\hat{\gamma}$ is slightly better.

Table V summarizes the lower bounds generated for problem 3 using the
procedure outlined in Section IV-B. Based on these results, the upper bounds
generated when $\gamma=\hat{\gamma}$ and $N$ = 2,500 (Table IV) are optimal at
99% confidence level. The upper bounds generated when $\gamma>\hat{\gamma}$
and $N$ = 400 (Table III) are 0.3 to 0.5% away from optimality at 99%
confidence level. Due to the high computational efforts to find feasible
solutions when $\gamma=\hat{\gamma}$ (with only marginal improvement to
solution quality), in the rest of our experiments we use $N=400$, $M=10$, and
$\hat{\gamma}$ = 0.05.

TABLE III: Upper Bounds: Problem 3,
$\gamma=0.10$, $\hat{\gamma}=0.05$

Sample | Process Time(hrs) | Solution Risk | No. of | Avg. Run
---|---|---|---|---
Size ($N$) | Min | Max | Avg | Min | Max | Avg | Feasible | Time (s)
50 | 8.97 | 9.00 | 8.99 | 0.000 | 0.132 | 0.021 | 7 | 739
100 | 8.95 | 8.98 | 8.97 | 0.000 | 0.141 | 0.016 | 9 | 1,010
150 | 8.95 | 8.98 | 8.97 | 0.000 | 0.123 | 0.016 | 6 | 1,059
200 | 8.93 | 8.98 | 8.95 | 0.000 | 0.118 | 0.016 | 7 | 1,208
250 | 8.93 | 8.98 | 8.96 | 0.000 | 0.120 | 0.014 | 8 | 1,241
300 | 8.97 | 8.98 | 8.97 | 0.000 | 0.116 | 0.012 | 9 | 1,353
400 | 8.93 | 8.98 | 8.96 | 0.000 | 0.097 | 0.013 | 10 | 1,646
600 | 8.93 | 8.97 | 8.96 | 0.000 | 0.097 | 0.013 | 10 | 2,176

TABLE IV: Upper Bounds: Problem 3, $\gamma=\hat{\gamma}=0.1$

Sample | Process Time(hrs) | Solution Risk | No. of | Avg. Run
---|---|---|---|---
Size ($N$) | Min | Max | Avg | Min | Max | Avg | Feasible | Time (s)
1000 | 8.93 | 8.93 | 8.93 | 0.000 | 0.105 | 0.013 | 7 | 2,744
2500 | 8.93 | 8.93 | 8.93 | 0.000 | 0.099 | 0.013 | 10 | 6,390
5000 | 8.93 | 8.93 | 8.93 | 0.000 | 0.100 | 0.013 | 9 | 12,715

TABLE V: Lower Bounds: Problem 3, $\gamma=\hat{\gamma}=0.1$

Sample | Lower | Gap(%) | Ave. Running
---|---|---|---
Size ($N$) | Bound(hrs) | Min | Max | Ave | Time (sec)
1000 | 8.93 | 0.00 | 0.00 | 0.00 | 2,744
2500 | 8.93 | 0.00 | 0.00 | 0.00 | 6,390
5000 | 8.93 | 0.00 | 0.00 | 0.00 | 12,715

##  VI Analysis of Results and Observations

Tables VI to IX summarize experimental results when solving problems 1 through
6 for different values of $q$, $\underline{q}$ and $\tau.$ Recall that
$\underline{q}$*100 and $q$*100 represent the minimum and maximum utilization
of the reactor. The value of $\tau$ represents the frequency by which the
chance constraint is satisfied. While constraint (26) enforces that the chance
constraint is met within the time horizon $\mathcal{T};$ constraints (34)
enforce that the chance constraint is met every $\tau$ periods. In our
experiments $\tau$ takes values 15 and 30 minutes.

$\displaystyle
Pr\left(\sum_{i\in\mathcal{I}_{r}}\sum_{s\in\mathcal{S}}\sum_{b\in\mathcal{B}}(f_{b}-f^{*})X_{isbt}\geq
0\right)$ $\displaystyle\geq 1-\gamma$ $\displaystyle\forall
t=\tau,2\tau,\ldots,\mathcal{T}.$ (34)

In our numerical analysis we report the value of the _coefficient of
variation,_ which measures the variability of biomass flow to the reactor, and
is calculated as the standard deviation over $t\in\mathcal{T}$ of:

$\frac{\sum_{i\in\mathcal{I}_{r}}\sum_{s\in\mathcal{S}}\sum_{b\in\mathcal{B}}X_{isbt}}{\sum_{i\in\mathcal{I}_{r}}\sum_{s\in\mathcal{S}}\sum_{b\in\mathcal{B}}\sum_{t\in\mathcal{T}}X_{isbt}}.$

Observation 1: _The processing rate of the reactor decreases, and the
processing time, the maximum inventory, and the coefficient of variation
increase as the requirements for meeting carbohydrate content become rigid._
The experiments summarized in Table VI considers $\tau=15$ and Table VII
considers $\tau=30.$ When $\tau=30$, the average processing time of problems 1
to 5 is 2.2% shorter, reactor’s processing rate is 2.5% higher, maximum
inventory level is 6.5% lower, coefficient of variation is 14.9% lower. A
similar observation is made when comparing the results of Tables VIII and IX.
When $\tau=30$, the average processing time of problems 1 to 5 is 1.6%
shorter, reactor’s processing rate is 1.6% higher, maximum inventory level is
1.4% lower, coefficient of variation is 4.5% lower. Furthermore, problem 6,
which represents the random sequence, could not find solutions which satisfied
the chance constraints for $\tau=15$ and $\tau=30$ minutes.

TABLE VI: Summary of Results: $\tau=15$ min, $\underline{q}=0.90,q=0.95$

Problem | Process Time | Reactor Rate | Flow to | Ave. Inv. | Max. Inv. | Coefficient
---|---|---|---|---|---|---
Nr. | (hrs) | (mg/h) | Reactor (mg) | (mg) | (mg) | Variation
1 | 13.62 | 2.26 | 30.75 | 1.12 | 4.54 | 0.0294
2 | 11.70 | 2.63 | 30.75 | 0.83 | 4.18 | 0.0454
3 | 8.96 | 3.43 | 30.75 | 0.15 | 0.87 | 0.0525
4 | 9.30 | 3.31 | 30.75 | 0.17 | 0.98 | 0.0223
5 | 13.25 | 2.32 | 30.75 | 0.98 | 4.22 | 0.0184
6 | Infeasible |  |  |  |  |

TABLE VII: Summary of Results: $\tau=30$ min, $\underline{q}=0.90,q=0.95$

Problem | Process Time | Reactor Rate | Flow to | Ave. Inv. | Max. Inv. | Coefficient
---|---|---|---|---|---|---
Nr. | (hrs) | (mg/h) | Reactor (mg) | (mg) | (mg) | Variation
1 | 13.62 | 2.26 | 30.75 | 1.12 | 4.54 | 0.0285
2 | 10.78 | 2.86 | 30.75 | 0.63 | 3.31 | 0.0460
3 | 8.96 | 3.43 | 30.75 | 0.14 | 0.83 | 0.0524
4 | 8.96 | 3.43 | 30.75 | 0.13 | 0.87 | 0.0524
5 | 13.25 | 2.32 | 30.75 | 0.98 | 4.28 | 0.0181
6 | 15.15 | 2.04 | 30.75 | 1.31 | 4.53 | 0.0239

TABLE VIII: Summary of Results: $\tau=15$ min, $\underline{q}=0.80,q=0.90$

Problem | Process Time | Reactor Rate | Flow to | Ave. Inv. | Max. Inv. | Coefficient
---|---|---|---|---|---|---
Nr. | (hrs) | (mg/h) | Reactor (mg) | (mg) | (mg) | Variation
1 | 13.45 | 2.29 | 30.75 | 1.12 | 4.54 | 0.0575
2 | 10.70 | 2.88 | 30.75 | 0.68 | 3.41 | 0.0965
3 | 8.93 | 3.44 | 30.75 | 0.19 | 1.23 | 0.1108
4 | 8.96 | 3.43 | 30.75 | 0.17 | 1.08 | 0.0955
5 | 11.82 | 2.60 | 30.75 | 0.74 | 3.20 | 0.0363
6 | Infeasible |  |  |  |  |

TABLE IX: Summary of Results: $\tau=30$ min, $\underline{q}=0.80,q=0.90$

Problem | Process Time | Reactor Rate | Flow to | Ave. Inv. | Max. Inv. | Coefficient
---|---|---|---|---|---|---
Nr. | (hrs) | (mg/h) | Reactor (mg) | (mg) | (mg) | Variation
1 | 13.42 | 2.29 | 30.75 | 1.11 | 4.54 | 0.0583
2 | 9.90 | 3.11 | 30.75 | 0.49 | 2.78 | 0.0985
3 | 8.93 | 3.44 | 30.75 | 0.19 | 1.22 | 0.1108
4 | 8.93 | 3.44 | 30.75 | 0.17 | 1.23 | 0.1107
5 | 11.82 | 2.60 | 30.75 | 0.74 | 3.49 | 0.0363
6 | 14.00 | 2.22 | 30.75 | 1.17 | 4.45 | 0.0524

Observation 2: _Deterministic modeling of the system leads to underestimation
of the processing rate of the reactor, processing time, inventory and
coefficient of variation._ Tables X and XI summarize the results from solving
the deterministic equivalent (the mean value problem) of problems 1 to 6. We
obtain the deterministic equivalent model by substituting constraint (26) with
constraints (35). In (35), $\bar{f}_{b}$ represents the average carbohydrate
rate of feedstock $b$. Since (35) are deterministic constraints, we do not
need to use the SAA method.

$\displaystyle\sum_{i\in\mathcal{I}_{r}}\sum_{s\in\mathcal{S}}\sum_{b\in\mathcal{B}}(\bar{f}_{b}-f^{*})X_{isbt}\geq
0\hskip 21.68121pt\forall t=\tau,2\tau,\ldots,\mathcal{T}.$ (35)

TABLE X: Summary of Deterministic Results: $\tau=15$ min,
$\underline{q}=0.90,q=0.95$

Problem | Process Time | Reactor Rate | Flow to | Avg Inv. | Max Inv | Coefficient
---|---|---|---|---|---|---
Nr. | (hrs) | (mg/h) | Reactor (mg) | (mg) | (mg) | Variation
1 | 9.00 | 3.42 | 30.75 | 0.14 | 1.09 | 0.0523
2 | 9.00 | 3.42 | 30.75 | 0.15 | 0.98 | 0.0522
3 | 8.93 | 3.44 | 30.75 | 0.10 | 0.72 | 0.0522
4 | 8.93 | 3.44 | 30.75 | 0.11 | 0.92 | 0.0525
5 | 13.25 | 2.32 | 30.75 | 0.98 | 4.54 | 0.0184
6 | 13.23 | 2.32 | 30.75 | 1.04 | 3.31 | 0.0118

TABLE XI: Summary of Deterministic Results: $\tau=30$ min,
$\underline{q}=0.90,q=0.95$

Problem | Process Time | Reactor Rate | Flow to | Avg Inv | Max Inv | Coefficient
---|---|---|---|---|---|---
Nr. | (hrs) | (mg/h) | Reactor (mg) | (mg) | (mg) | Variation
1 | 8.98 | 3.42 | 30.75 | 0.15 | 1.25 | 0.0525
2 | 8.98 | 3.42 | 30.75 | 0.15 | 1.05 | 0.0526
3 | 8.93 | 3.44 | 30.75 | 0.10 | 0.69 | 0.0523
4 | 8.93 | 3.44 | 30.75 | 0.11 | 0.91 | 0.0525
5 | 13.25 | 2.32 | 30.75 | 0.98 | 2.96 | 0.0183
6 | 13.23 | 2.32 | 30.75 | 1.05 | 4.52 | 0.0116

By comparing the results of Tables VI and VII with Tables X and XI,
respectively, we observe that the processing time of the deterministic model
is 9.67% shorter, reactor’s processing rate is 11.93% higher, maximum
inventory level is 89.28% lower and the coefficient of variation is 78.11%
lower then the corresponding results of the stochastic model.

Observation 3: _Short sequences that begin by processing bales with low and
medium moisture level perform best. These sequences lead to robust solutions._
Based on the results in Tables VI through IX, problems 3 and 4 perform the
best, and problem 3 outperforms problem 4. In both problems we use the short
feedstock-based sequence. For problem 3 the moisture-based sequence begins by
processing low moisture bales, followed by medium and high moisture level.
Such sequences allow the system to build up inventory which is used when
processing high moisture bales. Problems 1 and 2 perform worst then 3, 4 and
5. In these problems, we process every bale of a particular feedstock at a
time. Such a strategy leads to higher inventory level and longer processing
times. In many of the problem instances we solved, the Random sequence could
not find a feasible solution. In the cases when a feasible solution is found,
its quality is the worst than problems 1 to 5. This points to the importance
of sequencing bales to streamline biomass preprocessing.

Comparing the results of problem 3 in Tables VI through IX with Tables X and
XI one can see that the performance of problem 3 did not change. This
indicates that problem 3 is robust to variations in problem parameters.

In our numerical analysis we assume that the time unit $t=1$ minute
($t\in\mathcal{T}$). We resolve Problem 3 considering $t=0.25$ minutes. The
results are summarized in Table XII. Using smaller values of $t$ increase the
size of the problem and the running time of our proposed algorithms. However,
doing this leads to an 8.3% reduction of processing time of the system, 9%
increase of processing rate of the reactor, and 50% lower maximum inventory.

TABLE XII: Summary of Stochastic Results: $\tau=15$,
$\underline{q}=0.90,q=0.95$, $t=0.25$min

Problem | Process Time | Reactor Rate | Flow to | Ave. Inv. | Max. Inv. | Coefficient
---|---|---|---|---|---|---
Nr. | (hrs) | (mg/h) | Reactor (mg) | (mg) | (mg) | Variation
3 | 8.19 | 3.75 | 30.75 | 0.15 | 0.82 | 0.0524

Observation 4: _Increasing reactor’s reliability requires the use of high
quality biomass and additional inventory._ Tables VI and VII assume $q=0.95$
and $\underline{q}=0.90$, and Tables VIII and IX assume $q=0.90$ and
$\underline{q}=0.80$. The 5.56% increase of the average reactor utilization
($q$) from 0.90 to 0.95 led to 4.92% longer processing time, 2.44% decrease of
reactor’s processing rate, 6.64% increase of maximum inventory. We conclude
that a system can achieve a higher utilization rate of the reactor via
additional inventory level and high biomass quality.

Observation 5: _Chance constraint optimization modeling of the system
identifies robust and realistic solutions to our problem._ We resolve model
($Q$) as a stochastic optimization model. To accomplish this, we set
$\hat{\gamma}=0.00$, thus, constraint (26) becomes.

$\displaystyle
Pr\left(\sum_{i\in\mathcal{I}_{r}}\sum_{s\in\mathcal{S}}\sum_{b\in\mathcal{B}}(f_{b}-f^{*})X_{isbt}\geq
0\right)$ $\displaystyle\geq 1$ $\displaystyle\forall
t=\tau,2\tau,\ldots,\mathcal{T}.$ (36)

Since in (36) $P()\geq 1$, these constraints are to be satisfied all the time.
Thus, constraints (36) and (37) are equivalent. Thus, the solution of
($\hat{Q}$) for $\hat{\gamma}=0.00$ is the corresponding stochastic
optimization solution of ($\hat{Q}$).

$\displaystyle\sum_{i\in\mathcal{I}_{r}}\sum_{s\in\mathcal{S}}\sum_{b\in\mathcal{B}}(f_{b}-f^{*})X_{isbt}\geq
0$ $\displaystyle\forall t=\tau,2\tau,\ldots,\mathcal{T}.$ (37)

TABLE XIII: Summary of Stochastic Results: $\tau=15$ min,
$\underline{q}=0.90,q=0.95$

Problem | Process Time | Reactor Rate | Flow to | Avg Inv | Max Inv | Coefficient | Nr. of
---|---|---|---|---|---|---|---
Nr. | (hrs) | (mg/h) | Reactor (mg) | (mg) | (mg) | Variation | Feasible
1 | Infeasible |  |  |  |  |  | 0
2 | Infeasible |  |  |  |  |  | 0
3 | 10.21 | 3.05 | 30.75 | 0.47 | 2.08 | 0.0475 | 7
4 | 12.64 | 2.47 | 30.75 | 0.92 | 4.12 | 0.0276 | 7
5 | 13.41 | 2.30 | 30.75 | 1.01 | 4.30 | 0.0191 | 7
6 | Infeasible |  |  |  |  |  | 0

TABLE XIV: Summary of Stochastic Results: $\tau=30$ min,
$\underline{q}=0.90,q=0.95$

Problem | Process Time | Reactor Rate | Flow to | Avg Inv | Max Inv | Coefficient | Nr. of
---|---|---|---|---|---|---|---
Nr. | (hrs) | (mg/h) | Reactor (mg) | (mg) | (mg) | Variation | Feasible
1 | Infeasible |  |  |  |  |  | 0
2 | 14.17 | 2.17 | 30.75 | 1.25 | 4.54 | 0.0517 | 1
3 | 10.18 | 3.06 | 30.75 | 0.46 | 1.99 | 0.0477 | 7
4 | 10.43 | 2.99 | 30.75 | 0.48 | 2.15 | 0.0392 | 7
5 | 13.40 | 2.30 | 30.75 | 1.01 | 4.28 | 0.0199 7 |
6 | Infeasible |  |  |  |  |  | 0

Tables XIII and XIV summarize the results of the stochastic programming
formulation of ($\hat{Q}$). When $\tau=15,$ we did not find a feasible
solution in 3 out of 6 problems solved. Of the problems for which we found
feasible solutions, 7 out of the 10 replications were feasible. A similar
observation is made when $\tau=30.$ Furthermore, comparing the results of
Tables VI and VII with Tables XIII and XIV, we observe that the stochastic
model leads to 19.21% increase in processing time, 7.26% decrease in reactor’s
processing rate, 17.10% increase of maximum inventory, and 1.20% decrease in
the coefficient of variation.

##  VII Managerial Insights

###  VII-A Feasibility of Large-Scale Implementation of the Proposed Models

Successful large-scale implementation of the models developed depends on a
number of factors, such as, ($i$) standardizing the format of bales to
facilitate transportation and storage; ($ii$) use of Radio Frequency
Identification (RFID) technology to record and keep track of bale properties
from harvest to the biorefinery; ($iii$) use of sensing and real time
monitoring of material attributes, such as, moisture content since it changes
with time; ($iv$) use of automated material handling equipment to efficiently
move large quantities of biomass; ($v$) use of automated process control to
implement the sequencing algorithms we propose.

The focus of this research is on identifying biomass blends to meet process
requirements within a biorefinery, given the inventory of biomass bales. The
scope of this model could be extended to optimize the performance of the
overall supply chain. Such a model would focus on identifying how many bales
of different biomass feedstock to purchase given biomass availability in the
region, biomass price and quality. The results of this model can be used to
inform the design of long-term contracts among farmers and the biorefinery.

###  VII-B Simple Rules to Sequence Bales

The results of our numerical analysis point to the importance of sequencing
bales based on moisture level and carbohydrate content. Although in this paper
we present the results from solving 6 problems, each representing a different
sequence, we tested numerically a number of different sequences. We conducted
numerical tests, rather then solving an optimization model, for two main
reasons. First, integrating a sequencing model within our chance constraint
optimization would increase the computational challenge of solving the
problem. Second, since our model considers only 3 moisture levels and 4
feedstocks, we were able to identify numerically the characteristics of
sequences which work best. Below we provide simple rules that practitioners
can use to sort bales based on moisture level and carbohydrate contents.

Rule 1: _Sort Bales Based on Moisture Level_ : Let
$\mathcal{N}_{L},\mathcal{N}_{M},\mathcal{N}_{H}$ represent the number of
bales of low, medium and high moisture level. Let
$\psi(\mathcal{N}_{L},\mathcal{N}_{M},\mathcal{N}_{H})$ represent the greatest
common divisor (GCD) of $\mathcal{N}_{L},\mathcal{N}_{M},\mathcal{N}_{H}$.

_The rule is_ : Begin the sequence by processing
$\frac{\mathcal{N}_{L}}{\psi()}$ low moisture bales, followed by
$\frac{\mathcal{N}_{M}}{\psi()}$ medium moisture bales, followed by
$\frac{\mathcal{N}_{H}}{\psi()}$ high moisture bales. Repeat this sequence
$\psi()$ times to process every bale in the inventory.

In the case when the GCD of $\mathcal{N}_{L},\mathcal{N}_{M},\mathcal{N}_{H}$
does not exist, one can adjust these number by decreasing them by 1 or 2 bales
so that the resulting
$\mathcal{N}^{\prime}_{L},\mathcal{N}_{M}^{\prime},\mathcal{N}_{H}^{\prime}$
have a GCD. Next, the same Rule 1 applies for
$\mathcal{N}^{\prime}_{L}+\mathcal{N}_{M}^{\prime}+\mathcal{N}_{H}^{\prime}$
bales. The remaining bales are processed at the very end.

Example: Let $\mathcal{N}_{L}=24,\mathcal{N}_{M}=40,\mathcal{N}_{H}=16$ bales.
Thus, $\psi(24,40,16)=8.$ The sequence is $3L-5M-2H$, which repeats 8 times.

Rule 2: _Sort Bales Based on Biomass Characteristics_ : Let $\mathcal{N}_{q}$
represent the number of bales that meet carbohydrate requirements, and
$\mathcal{N}_{nq}$ represent the number of bales that do not. Let
$\psi(\mathcal{N}_{q},\mathcal{N}_{nq})$ be the GCD of $\mathcal{N}_{q}$ and
$\mathcal{N}_{nq}$.

_The rule is_ : Begin the sequence by processing
$\frac{\mathcal{N}_{q}}{\psi}$ bales, followed by
$\frac{\mathcal{N}_{nq}}{\psi}$ bales. Repeat this sequence $\psi$ times to
process every bale in the inventory.

In the case when the GCD of $\mathcal{N}_{q},\mathcal{N}_{nq}$ does not exist,
one can reduce $\mathcal{N}_{nq}$ to $\mathcal{N}^{\prime}_{nq}$ so that the
resulting $\mathcal{N}_{q},\mathcal{N}_{nq}^{\prime}$ have a GCD. Next, the
same Rule 2 applies for $\mathcal{N}_{q}+\mathcal{N}^{\prime}_{nq}$ bales. The
remaining ($\mathcal{N}_{nq}-\mathcal{N}_{nq}^{\prime}$) bales are processed
at the very end.

Example: Let $\mathcal{N}_{q}=20,\mathcal{N}_{nq}=60$ bales. Thus,
$\psi(20,60)=20.$ The sequence is $1q-3{nq}$, which repeats 20 times.

Rule 3: Let $f_{b}^{\gamma}$ represent the value of carbohydrate content which
corresponds to the $\gamma$-th percentile of biomass bales of feedstock $b$.
Let ${S}_{q}=\\{b\in\mathcal{B}|f_{b}^{\gamma}\geq f^{*}\\}$ and
${S}_{nq}=\\{b\in\mathcal{B}|f_{b}^{\gamma}<f^{*}\\}$. Let
$\underline{f}=\min_{b\in S_{q}}f_{b}^{\gamma}$, and let $N(\underline{f})$ be
the corresponding number of bales. Let $\overline{f}=\max_{b\in
S_{nq}}f_{b}^{\gamma}$, and let $N(\overline{f})$ be the corresponding number
of bales. The rule is outlined in Algorithm 2.

Algorithm 2 Rule 3 Procedure

Initialize: $S_{q}$, $S_{nq}$, $f^{*}$

1:Calculate: $\underline{f}$, $\overline{f}$, $N(\underline{f})$,
$N(\overline{f})$

2:Calculate: $\overline{d}=\underline{f}-f^{*}$ and
$\underline{d}=f^{*}-\overline{f}$;

3:While $\left(N(\underline{f})>0\And N(\overline{f})>0\right)$ Do

4: If  ($\overline{d}\geq\underline{d}$):

5: Sequence 2 bales of $S_{q}$ followed by
$\lfloor\frac{\overline{d}}{\underline{d}}\rfloor$ bales of $S_{nq}$

6: Let $N(\underline{f})=N(\underline{f})-2$,
$N(\overline{f})=N(\overline{f})-\lfloor\frac{\overline{d}}{\underline{d}}\rfloor$

7: Else: Sequence $\lceil\frac{\underline{d}}{\overline{d}}\rceil+1$ bales of
$S_{q}$ followed by 1 bale of $S_{nq}$

8: Let
$N(\underline{f})=N(\underline{f})-\lceil\frac{\underline{d}}{\overline{d}}\rceil-1$,
$N(\overline{f})=N(\overline{f})-1$

9:EndWhile

10:If $N(\underline{f})=0$

11: $S_{q}=S_{q}\setminus b$

12:Else

13: $S_{nq}=S_{nq}\setminus b$

14:If $\left(S_{q}\neq\emptyset\And S_{nq}\neq\emptyset\right)$ , GoTo line 1,
otherwise GoTo line 15

15:Add the remaining bales to the sequence; and return the sequence created.

Rule 4: _Sort Bales Based on Moisture Level and Biomass Characteristics_ :
This rule builds upon rules 1, 2 and 3.

_The rule is:_ ($i$) Sequence bales based on moisture level using rule 1.
($ii$) Sequence bales based on biomass characteristics using rules 2 or 3.
($iii$) Sequence bales in such a way that both sequences identified are
followed simultaneously.

Example: Let $\mathcal{N}_{L}=24,\mathcal{N}_{M}=40,\mathcal{N}_{H}=16$ bales
and $\mathcal{N}_{q}=20,\mathcal{N}_{nq}=60$ bales. Following rule 1 we get
$3L-5M-2H$. Following rule 2 we get $1q-3nq$. Let $S_{q}=\\{\text{Miscanthus,
Switchgrass}\\}$ and $S_{nq}=\\{\text{Corn stover 2-P, Corn stover 3-P.}\\}$

Figure 5: Sequence of problem 7.

Figure 5 demonstrates the sequence we created using the data in Section V-A.
This sequence is created using rule 4. We refer to this as problem 7. Notice
that, both sequences are followed. The shade of the table changes with the
moisture level to follow the pattern 3H-5M-2L. The values in the table change
to represent the feedstock used, and these feedstocks change based on 1q-3nq.
Table XV summarizes the results of solving problem 7. Comparing the
performance of problem 3 in Table VI with the results of Table XV, it is clear
that the performance of the sequence we develop using rule 4 is similar to the
performance of the short sequence (problem 3).

TABLE XV: Rule 4: Summary of Results: $\tau=15$ min,
$\underline{q}=0.90,q=0.95$

Problem | Process Time | Reactor Rate | Flow to | Ave. Inv. | Max. Inv. | Coef. of
---|---|---|---|---|---|---
Nr. | (hrs) | (mg/h) | Reactor (mg) | (mg) | (mg) | Var.
7 | 8.93 | 3.44 | 30.75 | 0.13 | 1.11 | 0.0514

##  VIII Summary and Conclusions

The objective of the proposed research is to develop process control
strategies to improve the performance of a biorefinery, and ensure a
continuous flow of biomass to the reactor, while meeting the requirements of
biochemical conversion process. The conversion process requires that the
carbohydrate content of the biomass fed to the reactor is at least 59.1%.
Satisfying this threshold is challenging because carbohydrate content is
stochastic. We propose a stochastic optimization model to identify the infeed
rate of biomass to the system, processing speed of equipment and inventory
level to ensure a continuous feeding of the reactor. We propose a sample
average approximation (SAA) method to approximate the model, and develop an
algorithm which solved the problem efficiently. The model is tested and
validated using data from the Process Development Unit of Idaho National
Laboratory.

We conduct extensive numerical analysis. The results of this analysis lead to
several observations, such as, the robustness of the short bale sequences we
propose. The corresponding process control leads to improvements on processing
time, processing rate of the reactor, inventory and coefficient of variations
(in feeding biomass to the reactor). Finally, the implementation of the models
developed requires the use of standardized bale format, Radio Frequency
Identification technology, sensing and real time monitoring of material
attributes, automated material handling equipment, and automated process
control. The scope of the model proposed can be extended to include the whole
supply chain. These models can be used to identifying how many bales of
different biomass feedstock to purchase given biomass availability in the
region, biomass price and quality. Thus, the results of this model can be used
to inform the design of long-term contracts among farmers and the biorefinery.

## References

  * [1] Biorefinery optimization workshop summary report, 2016.
  * [2] The renewable fuel standard (rfs): An overview, 2020.
  * [3] Warren P. Adams and Hanif D. Sherali.  Mixed-integer bilinear programming problems.  Mathematical Programming, 59(1):279–305, Mar 1993.
  * [4] Viknesh Andiappan, Michael Francis D Benjamin, Raymond R Tan, and Denny KS Ng.  Systematic framework for design of biomass energy systems: Addressing criticality via redundancy allocation.  Chemical Engineering Transactions, 61:787–792, 2017.
  * [5] Viknesh Andiappan, Michael Francis D Benjamin, Raymond R Tan, and Denny KS Ng.  An integrated framework to address criticality in biomass tri-generation systems via redundancy allocation.  Process Integration and Optimization for Sustainability, 3(1):65–73, 2019.
  * [6] Viknesh Andiappan and Denny KS Ng.  Synthesis of tri-generation systems: Technology selection, sizing and redundancy allocation based on operational strategy.  Computers & Chemical Engineering, 91:380–391, 2016.
  * [7] Timothy A Bell.  Challenges in the scale-up of particulate processes—an industrial perspective.  Powder Technology, 150(2):60–71, 2005.
  * [8] J.R. Birge and F. Louveaux.  Introduction to Stochastic Programming.  Springer, New York, 1997.
  * [9] Billal Nazim Chebouba, Mohamed Arezki Mellal, and Smail Adjerid.  Multi-objective system reliability optimization in a power plant.  In 2018 International Conference on Electrical Sciences and Technologies in Maghreb (CISTEM), pages 1–4, 2018.
  * [10] Chyi-Shin Chen, Akash Narani, Aigerim Daniyar, Joshua McCauley, Sarah Brown, Todd Pray, and Deepti Tanjore.  Ensemble models of feedstock blend ratios to minimize supply chain risk in bio-based manufacturing.  Biochemical Engineering Journal, page 107896, 2020.
  * [11] D.W. Coit and E. Zio.  The evolution of system reliability optimization.  Reliability Engineering and System Safety, 2018.
  * [12] Nathan C Crawford, Allison E Ray, Neal A Yancey, and Nick Nagle.  Evaluating the pelletization of “pure” and blended lignocellulosic biomass feedstocks.  Fuel Processing Technology, 140:46–56, 2015.
  * [13] Ryan Davis, Ling Tao, ECD Tan, MJ Biddy, GT Beckham, Christopher Scarlata, Jake Jacobson, K Cafferty, J Ross, J Lukas, et al.  Process design and economics for the conversion of lignocellulosic biomass to hydrocarbons: dilute-acid and enzymatic deconstruction of biomass to sugars and biological conversion of sugars to hydrocarbons.  Technical report, National Renewable Energy Lab.(NREL), Golden, CO (United States), 2013.
  * [14] D. Dentcheva, A. Prékopa, and A. Ruszczyński.  Concavity and efficient points for discrete distributions in stochastic programming.  Mathematical Programming, 89:55–77, 2000.
  * [15] S. D. Eksioglu, B. Gulcan, M. Roni, and S. Mason.  A stochastic biomass blending problem in decentralized supply chains.  Naval Research Logistics, 68(4):434–453, 2021.
  * [16] M. R. Garey and D.S. Johnsonand R. Sethi.  The complexity of flowshop and jobshop scheduling.  Mathematics of Operations Research, 1(2), 1976.
  * [17] Ronald L Graham, Eugene L Lawler, Jan Karel Lenstra, and AHG Rinnooy Kan.  Optimization and approximation in deterministic sequencing and scheduling: a survey.  In Annals of discrete mathematics, volume 5, pages 287–326. Elsevier, 1979.
  * [18] Berkay Gulcan, Sandra D Eksioglu, Yongjia Song, Mohammad Roni, and Qiushi Chen.  Optimization models for integrated biorefinery operations.  arXiv preprint arXiv:2101.03098, 2021.
  * [19] Lowell D Hill.  Grain grades and standards: Historical issues shaping the future.  University of Illinois Press, 1990.
  * [20] Haizhen Huang, Haibo Chen, Guohua Wang, and Jun Liu.  Combustion characteristics of coal and biomass blends and thermal dynamic analysis.  In 2009 International Conference on Energy and Environment Technology, volume 1, pages 555–559, 2009.
  * [21] Gurobi Optimization Inc.  Gurobi optimizer reference manual version 3.0.  2010\.
  * [22] J. Jacobson, M. Roni, P. Lamers, and K. Cafferty.  Biomass feedstock and conversion supply system design and analysis.  Technical report, Idaho National Lab.(INL), Idaho Falls, ID (United States), 2014.
  * [23] N. Jennett and M. Roni.  Inl dupont stover data report, 2020.
  * [24] Sujin Kim, Raghu Pasupathy, and Shane G Henderson.  A guide to sample average approximation.  Handbook of simulation optimization, pages 207–243, 2015.
  * [25] Fikri Kucuksayacigil, Mohammad Roni, Sandra D Eksioglu, and Qiushi Chen.  Optimal control to handle variations in biomass feedstock characteristics and reactor in-feed rate.  arXiv preprint arXiv:2103.05025, 2021.
  * [26] Kai Lan, Longwen Ou, Sunkyu Park, Stephen S Kelley, Burton C English, T Edward Yu, James Larson, and Yuan Yao.  Techno-economic analysis of decentralized preprocessing systems for fast pyrolysis biorefineries with blended feedstocks in the southeastern united states.  Renewable and Sustainable Energy Reviews, 143:110881, 2021.
  * [27] Kai Lan, Longwen Ou, Sunkyu Park, Stephen S Kelley, and Yuan Yao.  Life cycle analysis of decentralized preprocessing systems for fast pyrolysis biorefineries with blended feedstocks in the southeastern united states.  Energy Technology, 8(11):1900850, 2020.
  * [28] Kai Lan, Sunkyu Park, Stephen S Kelley, Burton C English, Tun-Hsiang E Yu, James Larson, and Yuan Yao.  Impacts of uncertain feedstock quality on the economic feasibility of fast pyrolysis biorefineries with blended feedstocks and decentralized preprocessing sites in the southeastern united states.  GCB Bioenergy, 12(11):1014–1029, 2020.
  * [29] M. Langholtz, B. Stokes, and L. Eaton.  2016 billion-ton report: Advancing domestic resources for a thriving bioeconomy, volume 1: Economic availability of feedstock.  2016\.
  * [30] J. Luedtke and S. Ahmed.  A sample approximation approach for optimization with probabilistic constraints.  SIAM Journal on Optimization, 19:674–699, 2008.
  * [31] Garth P McCormick.  Computability of global solutions to factorable nonconvex programs: Part i—convex underestimating problems.  Mathematical programming, 10(1):147–175, 1976.
  * [32] S. N. Mhasawade and L. A. Bewoor.  A survey of hybrid metaheuristics to minimize makespan of job shop scheduling problem.  In 2017 International Conference on Energy, Communication, Data Analytics and Soft Computing (ICECDS), pages 1957–1960, 2017.
  * [33] Akash Narani, NVSN Murthy Konda, Chyi-Shin Chen, Firehiwot Tachea, Phil Coffman, James Gardner, Chenlin Li, Allison E Ray, Damon S Hartley, Blake Simmons, et al.  Simultaneous application of predictive model and least cost formulation can substantially benefit biorefineries outside corn belt in united states: A case study in florida.  Bioresource technology, 271:218–227, 2019.
  * [34] A. Nemirovski and A. Shapiro.  Scenario approximations of chance constraints.  In G. Calafiore and F. Dabbene, editors, Probabilistic and Randomized Methods for Design under Uncertainty, pages 3–47. Springer, 1st edition, 2006.
  * [35] Guanqun Ni and Lei Chen.  Improved scheduling for the three-machine proportionate open shop and mixed shop minimum makespan problems.  IEEE Access, 8:186805–186812, 2020.
  * [36] Longwen Ou, Guanqun Luo, Allison Ray, Chenlin Li, Hongqiang Hu, Stephen Kelley, and Sunkyu Park.  Understanding the impacts of biomass blending on the uncertainty of hydrolyzed sugar yield from a stochastic perspective.  ACS Sustainable Chemistry & Engineering, 6(8):10851–10860, 2018\.
  * [37] V. Owens, D. Karlen, and J. Lacey.  Regional feedstock partnership summary report: Enabling the billion-ton vision.  2016\.
  * [38] David Padua.  Encyclopedia of parallel computing.  Springer Science & Business Media, 2011.
  * [39] Allison E Ray, Chenlin Li, Vicki S Thompson, Dayna L Daubaras, Nick Nagle, and Damon S Hartley.  Biomass blending and densification: impacts on feedstock supply and biochemical conversion performance.  Technical report, Idaho National Lab.(INL), Idaho Falls, ID (United States), 2017.
  * [40] M. Roni, D.N. Thompson, and D.S. Hartley.  Distributed biomass supply chain cost optimization to evaluate multiple feedstocks for a biorefinery.  Applied Energy, 254:113660, 2019.
  * [41] Mohammad S Roni, David Thompson, Damon Hartley, Erin Searcy, and Quang Nguyen.  Optimal blending management of biomass resources used for biochemical conversion.  Biofuels, Bioproducts and Biorefining, 12(4):624–648, 2018.
  * [42] M Sami, K Annamalai, and M Wooldridge.  Co-firing of coal and biomass fuel blends.  Progress in energy and combustion science, 27(2):171–214, 2001\.
  * [43] Bhavna Sharma, Craig Brandt, Devita McCullough-Amal, Matthew Langholtz, and Erin Webb.  Assessment of the feedstock supply for siting single-and multiple-feedstock biorefineries in the usa and identification of prevalent feedstocks.  Biofuels, Bioproducts and Biorefining, 14(3):578–593, 2020.
  * [44] Jian Shi, Kevin W George, Ning Sun, Wei He, Chenlin Li, Vitalie Stavila, Jay D Keasling, Blake A Simmons, Taek Soon Lee, and Seema Singh.  Impact of pretreatment technologies on saccharification and isopentenol fermentation of mixed lignocellulosic feedstocks.  BioEnergy Research, 8(3):1004–1013, 2015.
  * [45] Andrew S Tanenbaum and Maarten Van Steen.  Distributed systems: principles and paradigms.  Prentice-hall, 2007.
  * [46] Vicki S Thompson, John E Aston, Jeffrey A Lacey, and David N Thompson.  Optimizing biomass feedstock blends with respect to cost, supply, and quality for catalyzed and uncatalyzed fast pyrolysis applications.  BioEnergy Research, 10(3):811–823, 2017.
  * [47] Guido Van Rossum and Fred L Drake Jr.  Python tutorial.  Centrum voor Wiskunde en Informatica Amsterdam, The Netherlands, 1995\.
  * [48] Jingfu Wang, Lei Wen, Ying Chen, and Hengchao Han.  Study on the combustion characteristics and kinetics of biomass and coal char blended fuels.  In 2018 7th International Conference on Renewable Energy Research and Applications (ICRERA), pages 576–581, 2018.
  * [49] Wayne L Winston and Jeffrey B Goldberg.  Operations research: applications and algorithms, volume 3.  Thomson/Brooks/Cole Belmont^ eCalif Calif, 2004.
  * [50] Edward J Wolfrum, Nicholas J Nagle, Ryan M Ness, Darren J Peterson, Allison E Ray, and Daniel M Stevens.  The effect of biomass densification on structural sugar release and yield in biofuel feedstock and feedstock blends.  BioEnergy Research, 10(2):478–487, 2017.
  * [51] Jian Zhang, Guofu Ding, Yisheng Zou, Shengfeng Qin, and Jianlin Fu.  Review of job shop scheduling research and its new perspectives under industry 4.0.  Journal of Intelligent Manufacturing, 30(4):1809–1830, 2019.

## Acknowledgment

The funding of this work was provided by the U.S. Department of Energy, Office
of Energy Efficiency and Renewable Energy, Bioenergy Technologies Office under
award Number DE- EE0008255 and Department of Energy Idaho Operations Office
Contract No. DE-AC07-05ID14517. The authors greatly thank INL staff in the DOE
Biomass Feedstock National User Facility for providing technical data for this
study. Specifically, we acknowledge Neal A. Yancey and Jaya S. Tumuluru, who
are the individual INL contributors to this research who provided helpful
comments, data, and other forms of support for this analysis. The views
expressed herein are those of the author only, and do not necessarily
represent the views of DOE or the U.S. Government.

|  Dr. Dahui (Kevin) Liu received his Ph.D in industrial engineering from
Clemson University, Clemson, SC, USA, in 2020. Dr. Liu has worked for SAIC
Motor Corporation Limited, Shanghai Aircraft Manufacturing Corporation
Limited, and thyssenkrupp Uhde Chlorine Engineers (Shanghai) Corporation
Limited. He is currently a Postdoctoral fellow in the Industrial Engineering
Department, University of Arkansas. His research focuses in the area of
operation research with applications in transportation and energy.
---|---
|  Dr. Sandra D. Eksioglu is the John M. and Marie G. Hefley Professor in
Logistics and Entrepreneurship at the Department of Industrial Engineering,
University of Arkansas. Dr. Eksioglu’s expertise is in the areas of operations
research, network optimization, and algorithmic development with applications
in transportation, logistics, and supply chain.
---|---
|  Dr. Mohammad Roni –“Roni” to his friends and colleagues – is the Operations
Research Lead of Idaho National Laboratory’s Software and Data Sciences
department. In this position, he serves as principal investigator and
technical lead on research projects sponsored by the U.S. Department of
Energy. He is developing complex optimization models and algorithms applied to
supply chain network analysis, resiliency optimization, transportation and
infrastructure planning, and process control and design.
---|---

##  Appendix A Algorithms

Algorithm 3 Lower Bound Procedure

Initialize: $\delta=0.01$; $\hat{\gamma}>\gamma$;
$N\geq\frac{1}{2(\hat{\gamma}-\gamma)^{2}}\log{\big{(}\frac{1}{\delta}\big{)}}$;
$M=10$.

1:$\underline{\nu}\leftarrow\infty$

2:for $j=1,\ldots,M$ do

3: Sample $f_{bN}$ for $b\in\mathcal{B}$

4: Solve ($\hat{Q}$) using SAA Algorithm to find $\hat{\nu}_{Nj}$

5: if $\hat{\nu}_{Nj}\leq\underline{\nu}$ then

6: $\underline{\nu}\leftarrow\hat{\nu}_{Nj}$

7: end if

8:end for

9:Return $\underline{\nu}$

Algorithm 4 Upper Bound Procedure

Initialize: $\delta=0.01$; $C=0$; $\hat{\gamma}<\gamma$; $N$; d;
$N^{\prime}=10000$; $M=10$.

1:while $C<M$ do

2: $C\leftarrow 0$

3: for $j=1,\ldots,M$ do

4: Sample $f_{bN}$ for $b\in\mathcal{B}$

5: Solve ($\hat{Q}$) using SAA Algorithm to find ${\bf x}_{N}$

6: Sample $f_{bN^{\prime}}$ for $b\in\mathcal{B}$

7: Calculate $\bar{H}_{N^{\prime}}({\bf x}_{N},f)$

8: Calculate $\mathcal{U}_{N^{\prime}}({\bf x}_{N},f)$

9: if $\mathcal{U}_{N^{\prime}}({\bf x}_{N},f)\leq\gamma$ then

10: ${\bf x}_{N}$ is a feasible to ($Q$)

11: $C=C$ \+ 1

12: else

13: ${\bf x}_{N}$ is not feasible to ($Q$)

14: end if

15: end for

16: $N=N+d$

17:end while

18:Return $N,C$

Algorithm 5 Distributed and Parallel Computing Procedure

Initialize: $M=10$; $N^{c}\leq 8$; $P^{t}\leq 28$.

1:C $\leftarrow$ 1

2:for $j=1,\ldots,M$ do

3: if C $\leq N^{c}$ then

4: Assign $j$-th problem ($\hat{Q}$) to node C.

5: C $\leftarrow$ C \+ 1

6: else

7: C $\leftarrow$ 1

8: Assign $j$ problem ($\hat{Q}$) to node C.

9: C $\leftarrow$ C \+ 1

10: end if

11:end for

12:for $n=1,\ldots,N^{c}$ do

13: if Nr of problems assigned to n $>$ 1 then

14: Assign processors to problems

15: Solve ($\hat{Q}$) using Algorithm 1

16: else

17: Assign each of $P^{t}$ processors to the problem

18: Solve ($\hat{Q}$) using Algorithm 1

19: end if

20:end for

21:Return ${\bf x}_{N}$ for each of $M$ problems ($\hat{Q}$)

##  Appendix B Input Data

|  |  | Moisture Levels
---|---|---|---
|  |  | High | Medium | Low
Bale | Moisture (%) | 25 | 17.50 | 10
Dry bulk density (dry Mg/cubic meter) | 0.144 | 0.144 | 0.144
First-stage | Operating conditions | Screen size (mm) | 76.2 | 76.2 | 76.2
Dry bulk density (dry Mg/cubic meter) | 0.144 | 0.144 | 0.144
Moisture (%) | 25 | 17.5 | 10
Process performance | Energy consumption (kWh/dry Mg) | 25.68 | 21.95 | 7.10
Grinder | Moisture loss (%) | 4.77 | 3.00 | 0.5
Dry matter loss (%) | 1.50 | 1.50 | 1.50
Bulk density change (dry Mg/cubic meter) | -0.09 | -0.10 | -0.11
Maximum in-feed rate (dry Mg/hour) | 2.20 | 4.53 | 5.23
Separation | Operating conditions | Screen size (mm) | 6.35 | 6.35 | 6.35
Dry bulk density (dry Mg/cubic meter) | 0.053 | 0.041 | 0.039
Moisture (%) | 20.23 | 14.50 | 9.50
Unit | Process performance | Moisture loss (%) | 0.71 | 0.50 | 0.00
Bypass ratio (%) | 40.48 | 44.98 | 49.98
Second-stage | Operating conditions | Screen size (mm) | 11.11 | 11.11 | 11.11
Dry bulk density (dry Mg/cubic meter) | 0.053 | 0.041 | 0.039
Moisture (%) | 19.52 | 13.99 | 9.50
Process performance | Energy consumption (kWh/dry Mg) | 49.71 | 16.58 | 14.67
Grinder | Moisture loss (%) | 4.00 | 3.00 | 0.7
Dry matter loss (%) | 0.50 | 0.50 | 0.50
Bulk density change (dry Mg/cubic meter) | 0.066 | 0.082 | 0.090
Maximum in-feed rate (dry Mg/hour) | 1.59 | 2.80 | 5.23
Pelleting | Operating conditions | Dry bulk density (dry Mg/cubic meter) | 0.119 | 0.123 | 0.129
Moisture (%) | 15.52 | 10.99 | 8.80
Process performance | Energy consumption (kWh/dry Mg) | 90.39 | 60.63 | 55.12
Equipment | Moisture loss (%) | 3.92 | 1.50 | 0.00
Bulk density change (dry Mg/cubic meter) | 0.547 | 0.542 | 0.537
Maximum in-feed rate (dry Mg/hour) | 3.33 | 3.81 | 4.76
TABLE XVI: Technical process-related data for switchgrass TABLE XVII: Short Sequence Bale | Sequence Pattern
---|---
Nr. | Feedstock-Based | Moisture-Based
1-10 | (2S-1C2)-(2S-3C2)-(1S-1C2) | 3L-5M-2H
11-20 | (1M-2C2)-(2S-3C2)-(1S-1C2) | 3L-5M-2H
21-30 | (1S-2C2)-(1M-1S-3C2)-(1M-1C2) | 3L-5M-2H
31-40 | (2C2-1C3)-(3C2-2C3)-(1C2-1C3) | 3L-5M-2H
41-50 | (1M-2C3)-(3C2-2C3)-(1C2-1C3) | 3L-5M-2H
51-60 | (1M-2C2)-(2M-1C2-2C3)-(2C2) | 3L-5M-2H
61-70 | (2C2-1C3)-(2M-3C3)-(1M-1C3) | 3L-5M-2H
71-80 | (1C2-2C3)-(4C2-1C3)-(1C2-1C3) | 3L-5M-2H