Air Separation Unit

DEGREE PROJECT, IN CHEMICAL ENGINEERING , SECOND LEVEL
STOCKHOLM, SWEDEN 2015
Integrated Scheduling and Control of

an Air Separation Unit Subject to
Time-Varying Electricity Prices
TED JOHANSSON
KTH ROYAL INSTITUTE OF TECHNOLOGY
SCHOOL OF CHEMICAL SCIENCE AND ENGINEERING

Integrated Scheduling and Control of an Air Separation Unit
Subject to Time-Varying Electricity Prices
by
Ted Johansson
A Master’s Thesis at KTH Royal Institute of Technology
Supervisor: Dr. Michael Baldea*

Examiner: Dr. Per Alvfors**
July 13, 2015
*University of Texas at Austin, McKetta Department of Chemical Engineering
**KTH Royal Institute of Technology, Department of Chemical Engineering and Technology

Sammanfattning
I detta examensarbete presenteras en ny metod för att göra schemaläggningsbeslut inom processindustrin
och samtidigt ta hänsyn till det dynamiska beteendet hos processen. En modell av en luftseparerings-
anläggning som producerar kvävgas och utnyttjar ett rörligt elpris användes för att exemplifiera denna
metod. Modellen omfattade en kryogenisk destillationskolonn med en integrerad återloppskokare /kon-
densator, en multiströms värmeväxlare, en kompressor, två turbiner och en kondensator. Den innehöll
5079 ekvationer och 437 differentiella variabler. Dynamisk optimering användes för att approximera
det dynamiska beteendet hos processen vid skiftningar mellan olika driftpunkter. Den registrerade data
utnyttjades sedan för att identifiera en reducerad modell som fångade det transienta beteendet hos rele-
vanta processvariabler. Den reducerade modellen bestod av 525 ekvationer och 67 differentiella variabler.
Den identifierade modellen visade på god matchning mellan relevanta processvariabler i de simulerade
övergångarna och den reducerade modellen.
Den reducerade modellen användes för att optimera schemaläggningen av luftsepareringsprocessen så
att elkostnaden över en tredagars period minimerades. De optimala resultaten visade på en minskning av
kostnaden på 2.6 % jämfört med en konstant produktionstakt. Schemat implementerades och simulerades
i den fullt dynamiska modellen över de första 24 timmarna för att jämföra relevanta processvariabler med
den reducerade modellen. Resultaten visade på god matchning mellan de båda modellerna.
Detta examensarbete visar att en exakt reducerade modell kan användas för att snabbt hitta ett
optimalt schema över ett större processystem. Detta genom att kraftigt minska systemets storlek utan
att offra noggrannhet av det dynamiska beteendet.
Abstract
A novel framework for making plant scheduling decisions while considering the plant process dynamics
is presented in this thesis. A model of an air separation unit built to supply nitrogen gas and subject
to time-varying electricity prices was used to illustrate this framework. The model includes a cryogenic
distillation column with an integrated reboiler/condenser, a multi-stream heat exchanger, a compressor,
two turbines, and a liquefier. It consisted of 5079 equations and 437 differential variables. The dynamic
behavior of the process during operating point transitions was determined using dynamic optimization.
This data were used to establish a reduced order dynamic model of the system which captures the
transient behavior of relevant process variables. The reduced order model consisted of 525 equations and
67 differential variables. The identified model showed a good fit between the relevant process variables
in the simulated transitions and the reduced order model.
The air separation unit process schedule was optimized using the reduced order model to minimize
electricity cost over a three day time horizon. The optimal result showed a 2.6 % reduction in electricity
cost compared to a flat production rate. The optimal schedule was implemented and simulated in the
full dynamic model for the first 24 hours to compare the relevant process variables to the reduced model
predictions. The result displayed good match between the reduced model and the full dynamic model.
This thesis shows that an accurate reduced order dynamic model can be used for quickly finding
the optimal schedule of large process systems. This by greatly reducing the size and complexity of the
system without sacrificing accuracy of the dynamic behavior. Furthermore, it also shows the economic
benefits of the integrating scheduling and control to count for the dynamic behavior of the system.
Acknowledgments
I was fortunate enough to conduct my research for this thesis at the University of Texas at Austin in the
renowned McKetta Department of Chemical Engineering. I worked under the supervision and guidance
of Dr. Michael Baldea, who helped me achieve my goals in research and writing. This thesis would not
have been possible without his help. I would also like to extend my thanks and gratitude to his research
group, who were always welcoming and helpful. Richard Pattison and Cara Touretzky were especially
valuable resources, whose knowledge and experience were always readily provided to me. I am proud to
present this thesis, a product of a great team of researchers whose help and expertise must be recognized.
I would like to also acknowledge the strong foundation in chemical engineering that I built while at
KTH. Each of my professors imparted upon me the tools and knowledge I utilized in formulating and
carrying out my thesis work. Dr. Per Alvfors, as my examiner, was an invaluable link to my home
university and always available for counseling or advice. While I conducted my research abroad, it was
made possible by the amazing university that I can now call my alma mater.
Contents
1 Introduction 1
2 Background 3
2.1 Air Separation Unit . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3
2.1.1 Electricity Pricing Options for ASUs . . . . . . . . . . . . . . . . . . . . . . . . . . 3
2.1.2 Energy and Cost Reduction in ASUs . . . . . . . . . . . . . . . . . . . . . . . . . . 4
2.1.3 Dynamic Operation and Scheduling of ASUs Subject to Time-Varying Electricity
Prices . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4
2.1.4 Control Strategies of ASUs . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5
2.1.5 Dynamic models of ASUs . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5
3 Full Dynamic Model of an Air Separation Unit 7

3.1 Flow sheet overview . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7
3.2 Limitations and Assumptions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7
3.3 Model equations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8
3.3.1 Column . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8
3.3.2 Integrated Reboiler and Condenser . . . . . . . . . . . . . . . . . . . . . . . . . . . 11
3.3.3 Primary Heat Exchanger . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14
3.3.4 Compressor and Turbine . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21
3.3.5 Storage . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23
4 Integrated Scheduling and Control Using Reduced Order Models 26

4.1 Optimal Production Levels . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26
4.2 Dynamic Optimization to Find Optimal Transitions . . . . . . . . . . . . . . . . . . . . . 29
4.3 Reduced Order Model Identification . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30
4.4 Scheduling Using the Reduced Order Model . . . . . . . . . . . . . . . . . . . . . . . . . . 33
4.5 Implementation of Schedule on the Full Dynamic Model . . . . . . . . . . . . . . . . . . . 34
5 Results and Discussion 35

5.1 Model identification . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35
5.2 Scheduling Using the Reduced Order Model . . . . . . . . . . . . . . . . . . . . . . . . . . 38
5.3 Implementation of schedule in the full dynamic model . . . . . . . . . . . . . . . . . . . . 40
6 Conclusions and Future Work 44

6.1 Future Work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45
A Concepts 50
A.1 Differential Alegbraic Equation System . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50
A.2 Optimization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51
A.2.1 Optimization Problem Formulation . . . . . . . . . . . . . . . . . . . . . . . . . . . 51
A.2.2 Dynamic Optimization and Discretization Approaches . . . . . . . . . . . . . . . . 51
A.2.3 Solution Algorithms . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52
A.3 Reduced Order Modeling for Scheduling Applications . . . . . . . . . . . . . . . . . . . . . 52
A.3.1 Scheduling . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52
A.3.2 Reduced Order Modeling . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53
i
A.3.3 Model Identification . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53
B Results 55
B.1 Steady state optimization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55
B.2 Outputs From the Scheduling Using the Reduced Order Model . . . . . . . . . . . . . . . 55
ii
List of Figures
1.1 The hierarchy of decision making time scales in a chemical plant . . . . . . . . . . . . . . 2
2.1 Flow sheet of the process . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6
3.1 Overview of the column . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8

3.2 Overview of the integrated reboiler/condenser . . . . . . . . . . . . . . . . . . . . . . . . . 12
3.3 Sideview of the PHX with the turbine . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15
3.4 One segment in zone 1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16
3.5 Overview of the PHX . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17
3.6 Front view of the PHX . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18
3.7 A part of a layer in the PHX . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18
3.8 Zone 2 of PHX . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19
3.9 The liquefier and storage unit . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24
4.1 Scheduling procedure . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26

4.2 Block diagram of the transfer function. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32
4.3 Price schedule for April 11 - 13 2015 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33
5.1 Step changes in production set point considered during the dynamic optimization . . . . . 36
5.2 Data and transfer function models for step change responses. . . . . . . . . . . . . . . . . 37
5.3 Dew pressure data compared to transfer function model and non-linear model . . . . . . . 38
5.4 Schedule for production rate set point . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39
5.5 Hold up in storage unit . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39
5.6 Expected constraint and variable profiles from the dynamic scheduling using the reduced
order model for 72 hours . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40
5.7 Production of the plant using the dynamic and reduced order model for the first 24 hours
of the schedule found using the reduced order model. NRMSE between the reduced order
model and the full dynamic model is also given . . . . . . . . . . . . . . . . . . . . . . . . 41
5.8 Hold up in storage during 24 hours using the dynamic and reduced model . . . . . . . . . 41
5.9 Variables and constraints using the full dynamic model and the reduced order model. . . . 43
A.1 Examples of control profiles . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52

A.2 Transfer functions in series . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54
A.3 Transfer function in parallel . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54
A.4 Hammerstein-Wiener model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54
B.1 Optimal control variables for different production rates with polynomial approximations . 55
B.2 Expected variable and constraint outputs from the scheduling using the reduced order
model for 72 hours . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56
iii
List of Tables
3.1 Indices, variables and parameters used in the column model . . . . . . . . . . . . . . . . . 11

3.2 Indices, variables and parameters used in the integrated reboiler/condenser model . . . . 14
3.3 Indices, variables and parameters used in the PHX model . . . . . . . . . . . . . . . . . . 21
3.4 Indices, variables and parameters used in the compressor and turbine models . . . . . . . 23
3.5 Indices, variables and parameters used in the storage model . . . . . . . . . . . . . . . . . 25
4.1 Production levels used in the steady state optimization . . . . . . . . . . . . . . . . . . . . 27

4.2 Manipulated variables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30
4.3 Outputs in the scale bridging model and model order . . . . . . . . . . . . . . . . . . . . . 31
iv
Chapter 1
Introduction
Profitability of a process is dependent on minimizing the cost of raw materials and capital expenditures.
This is achieved through process integration, where material and energy streams are recycled to exploit
the efficient use of raw materials and capital equipment [1]. To further improve profitability, chemical
processes should be able to react to (potentially rapid) changes in market conditions by switching prod-
ucts, product purities, product grades, or production rates. These market conditions can include changes
in product demand and costs of raw materials (e.g., chemicals, utilities) [2].
A plant schedule identifies the order of products to be produced in a plant, including the purity
grades and production rates. The schedule is optimized to ensure that the product demands are met,
while minimizing the inventory cost [3]. On the other hand, process control is used to ensure that the
set points of the control variables (e.g. temperature, pressure, or concentration) are held in the presence
of disturbances (e.g. inlet purities, pressures, or temperatures) to minimize the production of products
that does not meet the specified requirements (off-spec products). Process control often uses a model
to describe the dynamic behavior of variables of interest in the process and ensure that operational
constraints are met. A control law is determined which gives a relationship between the manipulated
variables (e.g. valve position or heat removal rate) and the controlled variables. When a plant transitions
between operating states based on the schedule, the setpoints of the control variables are changed. The
controller adjusts the manipulated variables to track the setpoints and move the system to the new
steady state without violating process constraints [4, 5].
In the conventional chemical plant operating paradigm, scheduling and control are performed sepa-
rated over different time scales, see Figure 1.1. Scheduling decisions are typically made days in advance
and control moves take place over the course of minutes or hours [4]. When a chemical plant is rapidly
changing the product, product quality, or production rate to respond to market fluctuations, scheduling
decisions must be made in a similar time frame as the process dynamics. However, a static plant model,
which relies on a table of transitions, is typically used to determine the optimal operating schedule. It
is assumed that the lower level control layer can then guide these transitions exactly. The result is a
potentially suboptimal schedule because the dynamic behavior is not considered in the scheduling de-
cisions. The dynamic layers of the hierarchy should be considered to find the optimal schedule. This
requires the integration of the scheduling and control layers. One way to do this is to build a detailed
dynamic model of the system and then find the optimal transitions which minimize cost. Unfortunately,
this method is time consuming and requires extensive computational time [2, 3, 6, 7]. To overcome this
problem, a reduced order model can be used to make the optimal scheduling decisions as suggested by
Park et al. [8]. A reduced order model captures the dominant process dynamics accurately with much
fewer variables. Such models give an approximation of the process behavior and greatly reduce the
computational burden [2, 6].
A novel framework for making scheduling decisions using a reduced order dynamic model identified
from historical process data is established in this thesis. Here, an Air Separation Unit (ASU) case study
is considered to illustrate the benefits of this general framework. An ASU exhibits a multi-time scale
dynamic behavior during transitions between operating states and is a major electricity consumer [9,10].
It has been shown that the operating cost savings could be realized if an ASU is operated optimally
and change production rates in response to fluctuations in electricity price [11].In past works [11], the
savings were realized by reducing production rates and utilizing a storage unit to meet demand when
1
prices are high. Nevertheless, the optimal schedule was found using a static model. Without considering
the process dynamics, it is not possible to ascertain that the system is pushed to its transition limits, or,
conversely, that all of the suggested transitions are feasible.
Including accurate process dynamics in the scheduling decisions may improve the profitability of an
ASU. The operational cost is minimized by ensuring that the production rate changes are feasible and
that the most aggressive production rate changes are made. Determining the optimal schedule using a
full dynamic plant model often requires extensive computational power. It is also often too difficult to
solve over the (extended) time horizons of interest. However, if an accurate reduced order model of the
process can be identified with significantly fewer states, then the process dynamics could be accounted
for in the scheduling optimization, without increasing the computational burden.
In the literature, several works have applied reduced order models to optimal scheduling [2, 6, 8].
However, these works have investigated case studies of relatively small systems with few state variables.
An ASU is a much more complex system consisting of several interconnected unit operations. The ASU
model displays highly nonlinear behavior and has many state variables. In this work, a framework is
established for creating accurate data-driven models of large, complex process systems which capture the
dynamic behavior of the process variables of interest. In this framework, the models can be identified
from historical process data or by obtaining a dataset from an accurate dynamic plant model. A model
of an ASU was used to obtain a dataset which captures the dynamic behavior of the plant. The behavior
was described over a wide range of production rate changes such that an accurate reduced order model
could be established to be used for optimal scheduling.
Figure 1.1: The hierarchy of decision making time scales in a chemical plant [4]
Objectives
The goal of this thesis is to build an accurate dynamic model of an ASU with a liquid nitrogen storage
unit. The operational cost should then be minimized by operating the plant dynamically when subject
to time-varying electricity prices. A reduced order model is to be found which accurately represents the
dynamic behavior of process variables of interest. The predicted operational cost and dynamic behavior
of the two models should also be compared. The main objective is to evaluate the novel scheduling
approach that considers a simplified representation of the process dynamics.
This thesis is organized as follows: Chapter 1 and 2 gives a brief introduction and background.
Concepts used are presented in appendix A. This is followed by a detailed description of the mathematical
ASU model in Chapter 3. The full dynamic model of the ASU includes a cryogenic distillation column,
an integrated reboiler/condenser, a Primary Heat Exchanger (PHX), a compressor, a turbine, and a
liquefier. Chapter 4 describes the optimal scheduling methodology. The steady state plant operation
is first optimized at several production levels. Next, dynamic optimization is used to find the optimal
transitions between the different production levels. The optimal transitions are used to mimic historical
process data which in turn, are used to identify a reduced order model of the system that represents
the dynamic behavior of relevant variables in the plant. An optimal schedule is determined by dynamic
optimization using the reduced order model when subject to time-varying electricity prices. Finally, the
optimal schedule is implemented on the full dynamic model to determine the effectiveness of the identified
schedule. The results are presented and discussed in Chapter 5 with the conclusions and future work in
Chapter 6.
2
Chapter 2
Background
In this chapter, a brief overview of the concepts and theories used in the report are given. First, an account
over an air separation process is given together with the most used control and scheduling strategies.
Thereafter, a short summary of theoretical concepts regarding optimization and model identification is
provided.
2.1 Air Separation Unit

An ASU is used to separate the components of air and is used extensively in industry today to produce
high purity oxygen, nitrogen, and argon. Oxygen is mainly used in metal production and chemicals and
gasification. Gaseous nitrogen is used in the chemical and petroleum industries and argon is primar-
ily used for its inert properties [12]. For high product demands, cryogenic distillation is preferred to
adsorption methods for separating air. Cryogenic distillation consists of compressing air and cooling it
down to cryogenic temperatures (around 100 K). Separation of the components is then carried out in
a staged vapor-liquid equilibrium unit (a cryogenic distillation column). A multi-stream heat exchanger
and an integrated reboiler/condenser are used for heat integration [13]. Cryogenic air separation is an
energy intensive process and significant effort has been made to reduce the electricity demand required
for compression [10]. By as early as the 1970s, the first computer aided attempts to design and con-
trol ASUs were made. Slight design and operation changes to reduce energy consumption continue to
be implemented to this day. This has mainly been done by optimizing the heat and mass integration,
implementing more efficient machinery, and using advanced control strategies [12].
In cryogenic distillation, the transition times for many of the control variables are on the order of
hours. Adjustments in the production level must consider the dynamics of the process to ensure that
the transitions are feasible. Several authors have studied the dynamic behavior and optimal control of
ASUs during production rate transitions [12, 14–17].
2.1.1 Electricity Pricing Options for ASUs

An ASU has a very high electricity requirement to compress the incoming air. The American industrial
gas industry bought 19.4 TW h of electricity in 2010. This corresponds to about 2.5 % of all electricity
purchased by the manufacturing sector in the U.S that year [18]. Today, the most common electricity
pricing option for the industry is the time-of-use tariff where the price is set over a specified time horizon.
The price depends on the seasons of the year, the day of the week, and even the time of the day. The
price is set by the electricity producer based on a forecast of the generation cost and demand. This
option increases the risk for the producer, which results in higher prices paid by the customer. Another
pricing option is the real-time pricing tariff or time-varying electricity price. In this pricing scenario,
the price of electricity typically changes on an hourly basis which reflects current demand. While this
increases the risk for the customer, the base price is often lower than in the time-of-use pricing scenario.
This option could further help distribute the electricity usage in the grid over the duration of the day
(flattening the electricity demand curve) [10].
3
2.1.2 Energy and Cost Reduction in ASUs
Optimizing the design of an ASU aims to lower both the energy use of the plant and the overall cost.
Strategies to reduce the energy usage in the column include using heat pumps to increase thermal
efficiency [19] and structured packing to reduce pressure drop in the distillation column [20]. Aneke and
Wang [21] suggested using a Rankine cycle to generate electricity from the excess heat of the compressors.
Their results showed a possible 11 % reduction in the overall electricity requirement. Rizk et al. [22]
studied three different distillation column designs and compare their exergy losses. The results suggested
that the plant efficiency can vary significantly depending on the column design.
Research regarding the optimization of the design of ASUs have often utilized linear models [17, 19,
21–23]. However, the equations describing the ASU’s behavior are actually highly nonlinear. Detailed
nonlinear models should be used for accurate optimal design calculations. Zhu et al. [17] considered the
optimal design of ASUs under uncertainty with a detailed nonlinear model, and their results showed a
reduction in the capital cost of the plant.
The inlet air compressor is the main electricity consumer of an ASU. Improvements in the compres-
sor and turbine manufacturing process using computational fluid dynamics in combination with three
dimensional machining to find the optimal form of the different parts has increased the efficiency [24].
Turbines are used in an ASU to utilize pressure drops to generate electricity. Efficiency improvements
include the use of active magnetic bearings so that the frictional loss is reduced [25]. Multi-stream heat
exchangers are used to further reduce the energy consumption of the plant and allow for multi-stream
use [26]. Computer simulations further improve the design of the process by utilizing and integrating
the existing units [24].
ASUs are commonly integrated with their supplier and/or customers to reduce the energy consump-
tion and the capital cost [24]. Integrated gasification combined cycle, coal-based iron making, and
combined power or industrial gas production facilities can all benefit from integrating their usage of
oxygen or nitrogen with an ASU. The overall efficiency, power output, and cost saving can improve
compared to a non integrated process [27]. Elevated pressure (above 6.5 bar) ASUs could also be used,
which could lead to the usage of smaller and more compact cryogenic components and the reduction of
capital cost [13].
2.1.3 Dynamic Operation and Scheduling of ASUs Subject to Time-Varying

Electricity Prices
Several authors suggest that adjusting production rates to take advantage of time-varying electricity
prices can reduce the operating costs of an ASU due to the lower average prices [10, 11, 28, 29]. Optimal
process operation decisions must be made using demand and electricity price forecasts to minimize the
operating cost and increase the plant profits [30, 31]. Zhu et.al [30] investigated the optimal operation
during a 24-hour period with time-varying electricity prices and uncertain product demands. They
used a nonlinear static model with non-zero transitions approximated using linear relationships between
transition times and load changes. The 24-hour time horizon was divided into four periods with different
constant prices in each. They showed up to a 5 % decrease in operating costs during the 24-hour time
horizon. Furthermore, uncertainties in product demand could also be handled effectively with their
approach. Miller et al. [32] proposed a thermodynamic ideal work calculation to estimate the minimum
power requirement to find the maximum to minimum price ratios which determine when a plant changing
capacity to take advantage of time-varying electricity prices is profitable. Depending on the process and
economic assumptions, a ratio between two and seven (maximum to minimum electricity price) was
determined which made the scenario profitable.
Ierapetritou et al. [11] modeled an ASU with a liquid storage tank and assumed that the process was
subject to time-varying electricity prices. During periods of low electricity prices, the production rate
is increased beyond the product demand rate and the liquid nitrogen storage tank is filled. Conversely,
when the price is high, the production rate drops and the demand is fulfilled with inventory from the
storage tank. The results showed a substantial reduction in total electricity cost and encourage the
use of real time pricing. Furthermore, the authors suggested that this approach could help to flatten
the electricity demand curve throughout the day. Karwan and Keblis [10] used a linear process model
with discrete production levels to investigate when variable plant operation in response to time-varying
electricity price is profitable. Again, they use an optimization-based approach to find the schedule
4
that would minimize the operational cost. Their results showed that air separation plants with a great
flexibility in the production capacity benefit the most. However, time-varying electricity pricing has an
appeal even with a low degree of flexibility. Mitra et al. [31] investigated production scheduling under
time-varying demands and prices using a mixed-integer linear programming approach. They used a
deterministic surrogate model for the process based on the production mode, transition logic, and mass
balances. They defined each feasible production region by a convex region of already known feasible
operating points. Transition logics were used to handle the transition between regions. Their result
showed a 10 % improvement of the cost savings compared to a heuristic model. Storage capacity is also
shown to be an important factor in reducing operational cost. It is worth noting that most research to
date on ASU operation subjected to time-varying electricity prices has used simplified static linear plant
models to describe the system. These kinds of models do not capture the dynamic behavior of the plant
during production rate transitions, which remains an open research problem.
2.1.4 Control Strategies of ASUs

Pneumatic single loop controllers were mostly used to control ASUs before advanced computer systems
ever existed. These control systems can only operate on a single-input single-output base with necessary
rules, overrides, feedforward, etc. With more advanced computer systems, a shift towards using multi-
input, multi-output control strategies could be accomplished to minimize the operational cost. One such
strategy is the model predictive controller, which is based on an optimization algorithm to predict future
control moves of the manipulated variables to minimize or maximize a figure of merit. This method,
compared to single-input, single-output control system, has the advantage of improving the operating
performance. Model predictive control strategies are used extensively today in the air separation industry
and have been implemented by every major supplier of ASUs [12, 14].
Currently, many ASUs operate at a steady production rate with infrequent changes. Startups and
shutdowns are also uncommon. Linear dynamic models and model predictive controls adequately control
such systems to guarantee purity requirements [33]. However, when air separation plants shift pricing
options from time-of-use pricing to time-varying electricity pricing, the process nonlinearities become
clearer because of the frequent changes in production rates. This leads to increased requirements on
the control system to ensure purity demands. Nonlinear model-based control systems need to be put
in place to guarantee acceptable operation [34]. Several authors investigate the use of a non-linear
model predictive control system for an ASU [14, 33, 35, 36]. However, a typical first principles model
of a cryogenic distillation column is often too large and complex to solve on-line in a process. Zhu et.
al [34] investigated the use of a reduced order model representation of the distillation column. Their
results showed that the reduced model could predict composition responses adequately over a range of
different disturbances. Huang et. al [14], on the other hand, used an advanced nonlinear model predictive
controller that reduced the computational time of a cryogenic distillation column to one CPU second.
This sensitivity-based control strategy could also operate over a range of different operating conditions.
2.1.5 Dynamic models of ASUs

An ASU consists of several main components which affect the dynamic behavior and energy consumption
of the plant (see Figure 2.1). Components often included in a dynamic model are a cryogenic distillation
column, a PHX, a set of turbines and a compressor. A separate liquefier may also be included if a liquid
storage unit is considered in the model. The process dynamics are mainly influenced by the distillation
column and the PHX. The purity of the products in the column has response times on the order of hours
when subjected to production rate changes. This is due to the small temperature differences between
equilibrium stages in the column and the hold up in each tray [12]. The dynamic behavior is described
using mass and energy balances together with equilibrium relations and equations describing the physical
properties. More details are given in the next chapter.
5
N2 Gas Product
Integrated Reboiler/Condenser
Evaporator
Waste
Auxiliary Heat Exchanger
Storage
tank
Air Feed Liquefier Drain
Compressor
Storage Unit
Primary Heat Exchanger
Turbine 2
Expansion Valve
Column
Turbine 1
Figure 2.1: Flow sheet of the process
6
Chapter 3
Full Dynamic Model of an Air

Separation Unit
As described in the introduction, an ASU consists of the interconnection of several unit operations. Each
unit model is described by a system of Differential Algebraic Equations (DAE). In this work, gPROMS
4.0.0 [37] was used to construct the full dynamic model and formulate (and solve) the optimization
problems. In this chapter, an overview of the process flowsheet is described followed by a detailed
description of the model equations.
3.1 Flow sheet overview

Figure 2.1 shows the process diagram. The feed air is compressed from ambient conditions in a single
stage compressor. The compressed air is cooled to 300 K in a heat exchanger, before it enters the multi-
stream heat exchanger (PHX). The air is then cooled further in the PHX and a fraction of the cooled air
is withdrawn before the stream is liquefied. The gas fraction is expanded in a turbine and mixed with
the liquid fraction. The two phase air mixture is fed to the bottom of the cryogenic distillation column.
The column consists of 30 stages and an integrated reboiler/condenser. The pure nitrogen product at
the top of the column is split between the product stream and the reflux which is sent to the condenser
where the stream is liquefied and fed back to the top of the column. The liquid at the bottom of the
column is expanded adiabatically to further cool the stream and this stream is fed to the reboiler where
it is vaporized as it liquefies the mixture in the condenser. The waste vapor exits the reboiler and is sent
to the PHX along with the product stream. Adiabatic expansion of the product stream provides further
cooling, and the stream is re-passed through the PHX. Finally, excess product (beyond the gas demand)
is sent to a liquefier unit and stored as liquid nitrogen where it can be gasified to meet demand when
production rates drop.
3.2 Limitations and Assumptions

In this report, the mathematical models were limited to the main dynamic parts of the ASU: the PHX, the
column with an integrated reboiler/condenser, the compressor, the turbines and the liquefier. Parts not
included were the auxiliary heat exchangers, the pumps, the pipes, the valves, and the adsorption unit.
Furthermore, the column model was limited to only include the actual stages with an ideal efficiency,
and sumps were not included in the model. Other heat exchangers than the PHX were assumed to be
perfectly sized and use water as the cooling or heating source which was assumed to be free of charge.
The ASU was assumed to already exist with every equipment sized perfectly for a constant demand. The
capital cost was therefore not included in the calculations. Moreover, the plant was assumed to have a
constant demand profile from a single customer requiring nitrogen with a maximum allowable impurity
(oxygen and argon) of 1000 ppm. The air feed was assumed to consist of 78 % nitrogen, 21 % oxygen
and 1 % argon.
Other assumptions about the specific parts of the model are given in the text.
7
3.3 Model equations
In this section the equations describing the ASU model are presented. First, the dynamic models of
the column and the integrated reboiler/condenser are presented followed by the PHX model and the
compressor and turbine models. Lastly, the the liquefier and storage unit models are presented.
3.3.1 Column
Several authors have developed ASU column models [15,16,38]. In this work, the model has been derived
based on that developed by Huang et al. [14] with the following assumptions:
• Constant pressure drop on each stage
• Ideal vapor phases
• Well-mixed streams entering each stage
• Vapor liquid equilibrium is established on each tray
For every stage a DAE system was derived which included differential material and energy balances as
well as algebraic equations describing the equilibrium, hydraulics, and physical properties. An overview
of the column is shown in Figure 3.1.
Lcondenser , xcondenser,j
V1 , y1,j
T1 , P1 , M1
i=1
V2 , y2,j L1 , x1,j
Vi , yi,j Li-1 , xi-1,j
Ti , Pi , Mi i
Vi+1 , yi+1, j Li , xi,j
V30 , y30,j L29 , x29,j
F30 , xf30,j i = 30
T30 , P30 , M30
L30 , x30,j
Figure 3.1: Overview of the column
Overall mass balance

dMi
= Li−1 + Vi+1 − Li − Vi + Fi (3.1)
dt
8
where Mi is the molar hold up on stage i [mol], Li is liquid flow from stage i [mol s−1 ], Vi is vapor flow
from stage i [mol s−1 ] and, Fi is the feed flow at tray i [mol s−1 ] (which is nonzero only at the feed stage).
Component mass balance
dxi,j (Li−1 (xi−1,j − xi,j ) + Vi+1 (yi+1,j − xi,j ) − Vi (yi,j − xi,j ) + Fi (xfi,j − xi,j )
= (3.2)
dt Mi
where xi,j is the liquid mole fraction of component j at stage i, yi,j is the vapor mole fraction of com-
ponent j at stage i and xfi,j is the liquid mole fraction of component j in the feed.
Energy balance
f
dhL
i (Li−1 (hL L V L V L L
i−1 − hi ) + Vi+1 (hi+1 − hi ) − Vi (hi − hi ) + Fi (hi − hi )
= (3.3)
dt Mi
−1
where hLi is the molar enthalpy of the liquid at stage i [kJ mol ], hVi is the molar enthalpy of the vapor
f
at stage i [kJ mol ] and hi is the molar enthalpy of the feed stream [kJ mol−1 ]. The liquid and vapor
−1
enthalpies were computed as:

X
hL
i = (aL L
j Ti + bj )xj,i (3.4)
j
X
hVi = (aVj Ti + bVj )xj,i (3.5)
j
where aL L V V
j , bj , aj , and bj are constants for liquid and vapor specific to each component j and Ti is the
temperature at stage i [K]. aL L V V
j , bj , aj , and bj were calculated by data regressed from ASPEN [39]
Soave-Redlich-Kwong physical property method.
Summation equation X
1= yi,j (3.6)
j
Hydraulic equation
Li = kd Mi (3.7)
where kd is a tuning constant set to 0.5 [14].
Vapor-liquid equilibrium
yi,j pi = γi,j xi,j psat
i,j (3.8)
where pi is the total pressure on stage i [bar], psat
i,j is the saturation pressure of component j [bar] and
γi,j is the activity coefficient for component j. psat
i,j was calculated using Antoine’s equation:
βj
log psat
i,j = αj − (3.9)
δj + Ti
where αj , βj and δj are the Antoine’s coefficients for each component j [40] and Ti is the temperature at
each stage [◦C]. γi,j is the activity coefficient describing the nonideal behavior of the liquid for component
j. It was calculated using Margules equation:
AN2 ,O2 x2i,O2 + AN2 ,Ar x2i,Ar + (AN2 ,O2 + AN2 ,Ar − AO2 ,Ar )xi,O2 xi,Ar
γi,N2 = exp (3.10a)
RTi
AN2 ,O2 x2i,N2 + AO2 ,Ar x2i,Ar + (AN2 ,O2 + AO2 ,Ar − AN2 ,Ar )xi,N2 xi,Ar
γi,O2 = exp (3.10b)
RTi
AN2 ,Ar x2i,N2 + AO2 ,Ar x2i,O2 + (AN2 ,Ar + AO2 ,Ar − AN2 ,O2 )xi,N2 xi,O2
γi,Ar = exp (3.10c)
RTi
9
where Aj,k are the Margules constants for components j, k (j 6= k) [41] and R is the gas constant
[J mol−1 K−1 ].
Equations (3.1)-(3.10) described the dynamic model of a single stage i. However, the system was
an index 2 DAE system because the algebraic variable Vi could not be explicitly calculated from the
equations above.
The following procedure was used to reduce the index (notation is the same as presented by Huang
et.al [14]). As seen in equation (3.4) and (3.5), the enthalpy was a function of both temperature and
composition. Differentiating the left hand side of the energy balance, equation (3.3), gave [14]:
 
f
L X ∂hL Li−1 (hL L V L V L L
i−1 − hi ) + Vi+1 (hi+1 − hi ) − Vi (hi − hi ) + Fi (hi − hi )
 ∂hi T̄i + i
x̄i,j  = (3.11)
∂Ti j
xi,j Mi
dTi dxi,j ∂hL ∂hL

where T̄i = dt and x̄i,j = dt . ∂Ti
i
and i
xi,j can be explicitly calculated from equation (3.4) and (3.5):
∂hL
i
X
= xi,j aL
j (3.12a)
∂Ti j
∂hLi
= aL L
j Ti + bj (3.12b)
xi,j
An expression for T̄i could be found through derivation of the summation equation, equation (3.6), with
respect to time [14]:
X dyi,j X dKi,j dxi,j

0= = xi,j + Ki,j (3.13)
j
dt j
dt dt
γi,j psat
i,j
where Ki,j = pi is a function of temperature and composition. After applying the chain rule to
(3.13) and rearranging, an equation for dT dt could now be established [14]:
i
P h P ∂Ki,j i
dTi j x i,j x̄
k ∂xi,k i,k + K x̄
i,j i,j
T̄i = =− ∂Ki,j
(3.14)
dt
P
j xi,j ∂Ti
∂Ki,j ∂Ki,j
With ∂xi,k and ∂Ti being explicitly calculated as:
∂Ki,j psat
i,j ∂γi,j
= (3.15a)
∂xi,k pi ∂xi,k
∂Ki,j psat
i,j ∂γi,j γi,j ∂psat
i,j
= + (3.15b)
∂Ti pi ∂Ti pi ∂Ti
Now equation (3.11) is an algebraic expression containing Vi .
Vapor hold-up
Vapor holdup is usually assumed to be negligible [14] but to capture the actual process dynamics in
the model, the holdup was considered in this case. The vapor flow into every tray was calculated using
an artificial first order dynamic expression:
dVi 1
= (Vi+1 − Vi ) (3.16)
dt τ
where τ is a time constant equal to 30 s
Equation (3.1), (3.2), (3.4)-(3.12) and (3.14)-(3.16) were used to describe each stage i and component
j in the column model. 30 equilibrium stages were assumed in the column and they are numbered in a
10
increasing order starting from the top stage (1). The total pressure drop in the column was estimated
to 0.2 bar. The feed entered at the bottom of the column (stage 30). A nomenclature of the indices, the
variables, and the parameters used in the column model are given in Table 3.1
Table 3.1: Indices, variables and parameters used in the column model
Indices
i Tray number
j Component, j=1,2,3
k Component, k=1,2,3
L Liquid
V Vapor
f Feed
Variables
M Molar hold up
L Liquid flow
V Vapor flow
F Feed flow
x Molar liquid fraction
y Molar vapor fraction
γ Activity coefficient
h Molar enthalpy
p Pressure
psat Saturation pressure
T Temperature
dT
T̄ dt
dx
x̄ dt
Parameters Value Reference
a Enthalpy constant See reference [39]
b Enthalpy constant See reference [39]
kd Tuning constant 0.5 [14]
Aj, k Margulesconstant See reference [41]
α Antoine’s constant See reference [40]
β Antoine’s constant See reference [40]
δ Antoine’s constant See reference [40]
τ Time Constant 30 Estimated
R Gas Constant 8.314 [42]
3.3.2 Integrated Reboiler and Condenser

The model for the integrated reboiler and condenser was adapted from Cao [38]. It contained separate
models for the reboiler and condenser that are linked together via a heat transfer term. Figure 3.2
illustrates the reboiler/condenser model.
Condenser
The condenser model assumed zero material and energy holdup and that saturated liquid leaves the
condenser.
Mass Balance
0 = V1 rgasdraw − Vcondenser (3.17)
−1
where V1 [mol s ] is the top stage vapor rate, rgasdraw is the fraction of the top stage vapor flow that
is sent to the condenser and Vcondenser is the vapor molar flow to the condenser [mol s−1 ] which was
equivalent to the liquid molar flow leaving the condenser, Lcondenser
11
Vreboiler , yreboiler,j
Reboiler
Treboiler , Preboiler , Mreboiler
L30 , x30,j , preboiler
Qcondensor , UA
Condenser
Tcondenser , Pcondenser
Vcondenser , y1,j
Lcondenser , xcondenser,j
rgasdraw=Vcondenser/V1
Lreboiler , xreboiler,j Drain
Fproducts
V1 , y1,j L30 , x30,j , p30
Products, to PHX Waste, to PHX From top tray To top tray From bottom tray
Figure 3.2: Overview of the integrated reboiler/condenser
Energy Balance
0 = V1 rgasdraw hV1 − Lcondenser hL
condenser − Qcondenser (3.18)
−1
where hV1is the top tray molar vapor enthalpy [kJ mol ], hL
is the liquid molar enthalpy in the
condenser
condenser [kJ mol−1 ] and Qcondenser is the heat extracted from the condensation of vapor [kW]. The
liquid enthalpy was calculated using equation (3.4) and Qcondenser was calculated as:
Qcondenser = U A(Tcondenser − Treboiler ) (3.19)

−1
where U A is the product of the overall heat transfer coefficient and heat transfer area [kW K ], Tcondenser
is the temperature in the condenser [K] and Treboiler is the temperature in the reboiler [K].
Because all of the inlet vapor is condensed, the composition of the outlet liquid is the same as the
inlet vapor:
y1,j = xcondenser,j (3.20)
where y1,j is the top stage composition of component j and xcondenser,j is the liquid molar composition
in the condenser of component j.
Furthermore, because the liquid is assumed to be saturated, the condenser pressure can be calculated
using the saturation pressure of each component and the liquid composition:
X
pcondenser = xcondenser,j psat
j (3.21)
j
where pcondenser is the pressure of the condenser [bar] and psat

j is the saturation pressure of component
j [bar]. psat
j was calculated using equation (3.9).
Reboiler
The reboiler is modeled as an additional equilibrium stage with a heat input in the energy balance. The
expansion valve is modeled as a static adiabatic expansion, and was included in the reboiler section for
12
simplicity.
Mass Balance
dMreboiler
= L30 − Lreboiler − Vreboiler (3.22)
dt
where L30 is the liquid flow from the bottom tray in the column [mol s−1 ], Lreboiler is the drain flow
from the reboiler [mol s−1 ] and, Vreboiler is the waste vapor flow from the reboiler [mol s−1 ]. To keep the
molar hold up at a set point, a proportional controller with saturation was implemented on the drain
flow:
(
Lsp,reboiler − Kc (Msp,reboiler − Mreboiler ) if Mreboiler ≥ Msp,reboiler
Lreboiler = (3.23)
0 if Mreboiler < Msp,reboiler
where sp denotes set point and Kc is the controller gain.
Component Mass Balance
dxreboiler,j L30 (x30,j − xreboiler,j ) − Vreboiler (yreboiler,j − xreboiler,j )

= (3.24)
dt Mreboiler
where xreboiler,j is the liquid molar fraction of component j in the reboiler, x30,j is the liquid molar
fraction of component j in the bottom stage of the column and yreboiler,j is the vapor molar fraction of
component j in the reboiler.
Energy Balance
 
L L
 ∂hreb T̄reb +
X ∂hreb L30 (hL L V L
30 − hreb ) − Vreb (hreb − hreb ) + Qreb
x̄reb,j  = (3.25)
∂Treb j
xreb,j Mreb
where Qreb is the heat transferred from the condensor to the reboiler [kW] and is equal to Qcondenser ,
−1
hLreb is the liquid molar enthalpy in the reboiler [kJ mol ] calculated using equation (3.4) and hVreb is
−1
the vapor molar enthalpy in the reboiler [kJ mol ] calculated using equation (3.5), and hL 30 is the molar
enthalpy at the bottom of the column [kJ mol−1 ]. T̄reb and x̄reb,j were defined in an equivalent way as
for the column model, equation (3.11).
For the vapor liquid equilibrium, equation (3.8) was used in a similar way as for a tray to calculate
the composition of the waste stream. The pressure in the reboiler , preboiler was held constant at 2.5
bar. The indices, the variables and the parameters used in the integrated reboiler/condenser model are
given in Table 3.2.
13
Table 3.2: Indices, variables and parameters used in the integrated reboiler/condenser model
Indices
i Tray number
j Component, j=1,2,3
L Liquid
V Vapor
1 Top tray in column
30 Bottom tray in column
sp Set point
reb reboiler
Variables
M Molar hold up
L Liquid flow
V Vapor flow
x Molar liquid fraction
y Molar vapor fraction
γ Activity coefficient
h Molar enthalpy
p Pressure
T Temperature
UA Product of overall heat transfer coefficient and heat
transfer area
rgasdraw Fraction of product entering the condenser
Q Heat flow
dT
T̄ dt
dx
x̄ dt
Kc Controller gain 10−4 Estimated
preboiler Reboiler Pressure 2.5 [bar] Estimated
3.3.3 Primary Heat Exchanger

In this study, a brazed aluminum plate-fin multistream heat exchanger was assumed to be used as the
primary heat exchanger. It is built up of several layers of fins stacked together with a parting sheet
between them. A multistream heat exchanger is commonly used in plants where gas is separated into
its constituents, as an ASU. Its advantages are, among others, high surface area to volume ratio which
leads to smaller designs, ability to handle small temperature differences, and multistream use. [43]
The model of the PHX follows the work of Cao [38]. To simplify the model, the PHX is divided into
two zones. The first zone only considers the sensible heat change in the gas phase. The zone is further
divided into segments to accurately describe the temperature profile in the zone, see Figure 3.3. The air
feed is then divided into two streams. One stream is withdrawn from the heat exchanger as saturated
vapor and the other stream is liquefied in the second zone of the PHX. The temperature of air changing
phase in zone 2 is calculated as the average of the inlet and outlet temperature in the zone. The product
stream is fed to the column at high pressure and is withdrawn after the first zone and expanded and
re-passed through the heat exchanger as can be seen in Figure 3.3. The model of zones 1 and 2 consist
of differential equations describing the wall and stream temperature and algebraic equations describing
the physical properties of the PHX and the streams.
14
Turbine
High
pressure
Low product
Pressure
product
Segment 1 Segment n Segment 50 Segment 1

Waste
Air Feed Liquid air feed
Vapor air feed
Zone 1 Zone 2
Figure 3.3: Sideview of the PHX with the turbine
Zone 1
Zone 1 consisted of differential equations describing the wall and stream temperature in every segment.
In each segment, the streams were assumed to be well mixed and the wall temperature uniform. Figure
3.4 shows an overview of a segment in zone 1.
15
Segment n
F1,1 , hV1,1,n F1,1 , hV1,1,n+1

E1,1,n , T1,1,n , P1,1,n , UA1,1,n , V1,1,n , Tm1,n
F1,2 , hV1,2,n F1,2 , hV1,2,n+1

E1,2,n , T1,2,n , P1,2,n , UA1,2,n , V1,2,n , Tm1,n
F1,3 , hV1,3,n F1,3 , hV1,3,n+1

E1,3,n , T1,3,n , P1,3,n , UA1,3,n , V1,3,n , Tm1,n
F1,4 , hV1,4,n-1 F1,4 , hV1,4,n

E1,4,n , T1,4,n , P1,4,n , UA1,4,n , V1,4,n , Tm1,n
Figure 3.4: One segment in zone 1
Metal wall temperature

m
dTk,n X
mk,n Cpmetal = m
U Ak,j,n (Tk,j,n − Tk,n ) (3.26)
dt j
where mk,n is the metal wall mass in zone k and segment n [kg], Cpmetal is the specific heat capacity
of the metal wall [J kg−1 K−1 ], Tk,j,n is the temperature in zone k and segment n of stream j [K], Tk,nm
is the metal wall temperature in zone k and segment n, and U Ak,j,n is the product of the overall heat
transfer coefficient and the heat transfer area in zone k and segment n for stream j [J K−1 ]. The streams
were numbered in increasing order as: the high pressure product stream, low pressure product stream,
the waste stream, and the feed air stream. The metal wall mass was calculated as:
rk mmetal
mk,n = (3.27)
Nsegments
where rk is the fraction of the PHX that is in zone 1, Nsegments is the number of segments in zone 1 and
mmetal is the total metal mass in the PHX [kg] calculated as:
NS
!
X
mmetal = ρmetal W LH − Vj (3.28)
r
where ρmetal is the metal density [kg m−3 ], W, L and H are the width, length and height of the PHX [m]
(see Figure 3.5), and Vj is the total volume occupied by stream j in the PHX [m3 ]. Vj was estimated
from the volume occupied by stream j in each slot:

1 1
Vjslot = bj − bj tj + − tj tj (L − 2tedge ) (3.29)
Nj Nj
16
H
Figure 3.5: Overview of the PHX
where Vjslot is the volume occupied by stream j in each slot [m3 ], N1j is the width of a single slot for
stream j [m], bj is the distance between the parting sheets for stream j [m], tj is the fin thickness for
stream j [m] and tedge is the thickness of the PHX wall [m] (see Figure 3.6 and 3.7).
The volume occupied by stream j in each layer was then calculated as:
Vjlayer = (W − 2tedge )Nj Vjslot (3.30)
Finally, the total volume occupied by each stream in the PHX could then be estimated as:
Vj = Vjlayer nj (3.31)
where nj is the number of layers for each stream.
Energy Balance
The energy balance in each segment of zone 1 was:
dE1,j,n
= F1,j (hV1,j,n+1 − hV1,j,n ) + U A1,j,n (T1,n
m
− T1,j,n ) (3.32)
dt
where E1,j,n is the energy hold up in stream j in segment n and zone 1 [kJ], Fk,j is the molar flow of
stream j in zone 1 [mol s−1 ], and hV1,j,n is the vapor enthalpy of stream j in segment n [kJ mol−1 ]. The
enthalpy of the streams in each segment was calculated using equation (3.5). E1,j,n was defined as:
E1,j,n = hV1,j,n V1,j,n ρV1,j,n (3.33)

3
where V1,j,n is the volume occupied in each segment by stream j [m ] and was calculated from:
r1 Vj
V1,j,n = (3.34)
Nsegment
The molar density of stream j in zone 1 and segment n, ρV1,j,n [mol m−3 ], in equation (3.33) was calculated
from the ideal gas law:
p1,j,n
ρV1,j,n = (3.35)
RT1,j,n
17
Figure 3.6: Front view of the PHX
Figure 3.7: A part of a layer in the PHX
where p1,j,n is the pressure of stream j in segment n and zone 1 [Pa], T1,j,n is the temperature of stream
j in zone 1 and segment n [K], and R is the gas constant [J mol−1 K−1 ]. Zone 1 was divided into 50
segments.
Zone 2
The air feed stream, stream 4, was divided at the end of zone 1. One fraction was extracted as vapor
and the other continued into zone 2 and phased changed into liquid air, see Figure 3.8. To simplify the
calculations only one segment was considered in zone 2 i.e Nj,segment = 1. The temperature for each
stream in zone 2 was calculated as an average of the inlet and outlet temperature of that stream:
18
Segment 1
F1,1 , hV2,1,out , T1,out F2,1 , hV2,1,in , T1,in

E2,1 , T2,1 , P2,1 , UA2,1 , V2,1 , Tm2

E2,2 , T2,2 , P2,2 , UA2,2 , V2,2 , Tm2

E2,3 , T2,3 , P2,3 , UA2,3 , V2,3 , Tm2
F1,4 F2,4 , hV2,4,in , T4,in F2,4 , hL2,4,out , T4,out

E2,4 , T2,4 , P2,4 , UA2,4 , V2,4 , Tm2
kwithdrawn=F2,4/F1,4
Fvapor
Zone 1 Zone 2
Figure 3.8: Zone 2 of PHX
Tj,in + Tj,out
T2,j = (3.36)
2
where Tj,in is the inlet temperature of stream j and Tj,out is the outlet temperature of stream j.
The change in metal wall temperature was captured by the same equation as in zone 1, equation
(3.26). The energy balance was formulated equivalent to that of zone 1 for streams 1-3:
dE2,j
= F2,j (hV2,j,in − hV2,j,out ) + U A2,j (T2m − T2,j ) (3.37)
dt
where hV2,j,in is the vapor enthalpy of stream j at the inlet of zone 2 [kJ mol−1 ] and hV2,j,out is the vapor
enthalpy of stream j at the outlet of zone 2 [kJ mol−1 ]. The inlet and outlet enthalpies were calculated
using equation (3.5) and the respective inlet and outlet temperature.
A fraction of stream 4 changed phase from gaseous air into liquid air and the energy balance was
calculated as:
dE2,4
= F2,4 (hV2,4,in − hL m
2,4,out ) + U A2,4 (T2 − T2,4 ) (3.38)
dt
where hliq
2,3,out is the liquid enthalpy of stream 4 at the outlet of zone 2 [kJ mol
−1
] calculated using
equation (3.4) and the outlet temperature of that stream. F2,4 is the liquid air feed stream [mol s−1 ]
calculated as:
F2,4 = kwithdrawn F1,4 (3.39)

where kwithdrawn is the fraction of the feed flow that was withdrawn from the PHX as vapor prior to
zone 2 and F1,4 is the molar air stream in zone 1 [mol s−1 ].
E2,j was defined in the same way as in equation (3.33) and the average temperature in zone 2, T2,j .
E2,4 was calculated as:
19
E2,4 = hL L
2,4 V2,4 ρ2,4 (3.40)
−1
where hL2,4 is the liquid enthalpy for stream 4 in zone 2 [kJ mol ] calculated using equation (3.4) and
the average temperature, T2,4 . V2,4 is the volume occupied by stream 4 in zone 2 [m3 ] calculated in the
−3
same way as in zone 1, equation (3.34). ρL 2,4 is the liquid molar density of stream 4 in zone 2 [mol m ]
and was calculated using the Rackett equation [38, 44]:
" 2/7 #
RTic

1 T2,4
ln = ln c + 1+ 1+ c ln(ZiRa ) (3.41)
ρLi p i T i
NC
X
ρL
2,4 = xi ρL
i (3.42)
i=1
−3
where ρLi is the liquid molar density of component i in stream 4 [mol m ], Tic is the critical temperature
for component i in stream 4 [K] [45], pi is the critical pressure for component i in stream 4 [Pa] [45], ZiRa
c
is the Rackett compressibility factor of component i in stream 4 [46], and xi is the liquid molar fraction
of component i in stream 4.
mk,n and ρV2,j was calculated in the same way as in zone 1 but with only 1 segment (equation (3.27)
and (3.35), respectively). Variables, indices and parameters used in the PHX model are given in Table
3.3.
20
Table 3.3: Indices, variables and parameters used in the PHX model
Indices
k Zone
n Segment
j Stream
i Component
L Liquid
V Vapor
in Inlet
out Outlet
Variables
Tm Metal wall temperature
T Stream temperature
mk,n Metal wall mass in zone k, segment n
mmetal Total metal mass
Vj Volume occupied by stream j
Vslot
j Volume occupied by stream j in each slot
Vlayer
j Volume occupied by stream j in each layer
E Energy hold up
F Molar flow
h Molar enthalpy
ρ Molar density
p Pressure
kwithdrawn Fraction liquid air entering zone 2
x Molar fraction
T Temperature
Cp Specific heat capacity of metal wall 892 [J kg−1 K] Estimated
UA1,j,n Product of overall heat transfer coefficient and heat 40000 [J K−1 ] Estimated
transfer area in zone 1
UA2,j,n Product of overall heat transfer coefficient and heat 10000 [J K−1 ] Estimated
transfer area in zone 2
rk Fraction of PHX that is in zone 1 0.954817 Estimated
Nk,segments Number of segments in zone 1 50 Estimated
ρmetal Metal density [kg m−3 ] 2731 Estimated
W Width of PHX [m] 1.5 Estimated
L Length of PHX [m] 8.2 Estimated
H Height of PHX [m] 3 Estimated
Nj Number of fins of stream j [70 70 70 100] Estimated
bj Distance between parting sheets for stream j [m] 6.75e-3 Estimated
tj Fin thickness for stream j [m] 0.4e-3 Estimated
tedge PHX wall thickness [m] 1e-2 Estimated
nj Number of layers for each stream 59 Estimated
R Ideal gas constant 8.314 [42]
Tci Critical temperature of component i See reference [45]
pci Critical pressure of component i See reference [45]
ZRa
i Rackett compressibility factor of component i See reference [46]
3.3.4 Compressor and Turbine

The model of the compressor follows the developments in Green and Perry [44]. The calculations were
based on a polytropic compression and the power demand from the compressor was calculated as:
21
GHpoly
W = (3.43)
ηpoly
where W is the power demand [W], G is the mass flow [kg s−1 ], Hpoly is the polytropic head [J kg−1 ] and
ηpoly is the polytropic efficiency. Hpoly was calculated as:
" (n−1)/n #
n pout
Hpoly = ZRTin −1 (3.44)
n−1 pin
where n is the polytropic exponent, Z is the compressibility factor, R is the gas constant [J mol−1 K−1 ],
Tin is the temperature of the inlet stream [K], and pout and pin are the pressure exiting and entering the
compressor [bar], respectively. ηpoly in equation (3.43) was calculated as:
(k − 1)/k
ηpoly = (3.45)
(n − 1)/n
where k is defined as the ratio of specific heats:
Cp
k= (3.46)
Cv
Values for Cp and Cv was found in CRC handbook of chemistry and physics [42] for oxygen and nitrogen
and are assumed to be independent of temperature.
When a fluid is compressed the relationship between pressure, p, and Volume, V , can be described
using the polytropic exponent:
pV n = Constant (3.47)
For a polytropic process the exponent is constant and equation (3.47) can be rearranged as follows:
ln p + n ln V = Constant (3.48)
An empirical relationship can be used if no data are available. Perry and Green [44] suggest a relationship
between the adiabatic efficiency and the polytropic exponent:
pout (k−1)/k
pin −1
ηad = pout (n−1)/n
(3.49)
pin −1
where ηad is the adiabatic efficiency. It was approximated to 80 % in this study.
Turbine
The turbine model was based on the same equations as for the compressor and the model was used for both
turbines in the process flowsheet. It was assumed that the polytropic exponent could be approximated
in the same way as in equation (3.49) and that the adiabatic efficiency was the same for the turbine as
for the compressor. The power supply from the turbine was then calculated as:
W = GHpoly ηpoly (3.50)

Indices, variables and parameters used in the compressor and turbine models are presented in Table 3.4
22
Table 3.4: Indices, variables and parameters used in the compressor and turbine models
Indices
poly Polytropic
in Inlet
out Outlet
ad Adiabatic
Variables
W Power
G Mass flow
H Polytropic head
η Efficiency
p Pressure
T Temperature
k Ratio of specific heats
n Polytropic exponent
Z Compressibility factor 1 Estimated
R Ideal gas constant 8.314 [42]
Cp,air Specific heat capacity for air [kJ kg−1 K] 1.007 [42]
Cp,N2 Specific heat capacity for nitrogen [kJ kg−1 K] 1.123 [42]
Cv,air Specific heat capacity for air [kJ kg−1 K] 0.7181 [42]
Cv,N2 Specific heat capacity for nitrogen [kJ kg−1 K] 0.7710 [42]
ηad Adiabatic efficiency [%] 80 Estimated
3.3.5 Storage
The storage unit consisted of a liquefier, a storage tank, and an evaporator. The liquefier was approx-
imated by a steady state model as the dynamics are typically faster than the ASU. Furthermore, the
evaporator was assumed to work with an external water cooled heat exchanger and was not included as
a model, see Figure 3.9.
Liquefier
The liquefier was modeled as an ideal cycle with an isothermal compressor and an isentropic turbine at
steady state [47], see Figure 3.9. The energy balance was derived as:
Wc + F h1 = WT + Qr + F hf (3.51)
where Wc is the work required by the compressor [W], WT is the work done by the turbine [W], Qr is
the heat produced by the compression [W], F is the molar flow to the liquefier [mol s−1 ], h is the molar
enthalpy [J mol−1 ]. Rearranging equation (3.51) gave the total energy consumption of the liquefier:
Wnet = Wc − WT = Qr − F (h1 − hf ) (3.52)

The liquid at the outlet of the liquefier is assumed to be saturated and at the same pressure as the inlet.
The second law of thermodynamics was used to calculate Qr . An entropy balance around the system
was established:
Qr
F s1 = F sf + (3.53)
T1
Rearranging gives:
Qr = T1 F (s1 − sf ) (3.54)
−1 −1 −1 −1
where s1 is the entropy at the inlet [J mol K ], sf is the entropy at the outlet [J mol K ], and T1
is the temperature at the inlet [K].
23
Gaseous N2 in
F, h1 , T1 , s1
Gaseous N2 out
Evaporator
Compressor
Turbine Fout
F, hf , sf
Liquid N2 Storage tank
Liquefier
Figure 3.9: The liquefier and storage unit
Combining equation (3.51) and (3.53) gave an equation for the overall work required by the liquefier:
Wnet = T1 F (s1 − sf ) − F (h1 − hf ) (3.55)

The real work required by the liquefier was given by:
Wnet
Wreal = (3.56)
η
where η is the overall efficiency of the system and was estimated to be 40 %.
The amount of liquid sent to or withdrawn from the liquefier was based on the product demand:
(
Fproduced − Fdemand if Fproduced ≥ Fdemand
F = (3.57)
Fdemand − Fproduced if Fproduced < Fdemand
where Fproduced is the molar flow of nitrogen produced from the ASU [mol s−1 ] and Fdemand is the product
demand [mol s−1 ]. The demand was assumed to be constant at 20 mol s−1 corresponding to about 50
t day−1 .
The flow of nitrogen produced changed in response to the variability of the electricity prices. To
ensure a constant flow of nitrogen product out from the system, the extra nitrogen produced at low
prices was sent to the liquefier and stored as liquid product in the storage tank. When prices are high,
stored liquid nitrogen was evaporated to satisfy the demand.
24
Storage Tank
The liquefied nitrogen was stored in a cryogenic storage tank at the same pressure as the inlet pressure
to the liquefier. The tank was assumed to be well insulated and no heat was neither added nor lost from
the tank. It was modeled as an accumulation term:
dMstorage
= F − Fout (3.58)
dt
where Mstorage is the molar holdup in the tank [mol] and Fout is the molar flow out from the tank
[mol s−1 ]. The liquid flow out from the tank was evaporated in an external heat exchanger to the same
temperature as the product flow from the primary heat exchanger.
The indices, variables, and parameters used in the storage tank model are presented in Table 3.5
Table 3.5: Indices, variables and parameters used in the storage model
Indices
c Compressor
T Turbine
f Outlet
1 Inlet
net Net work
Variables
W Power
Qr Heat from compression
h Enthalpy
F Molar flow
s Entropy
k Ratio of specific heats
n Polytropic exponent
Mstorage Hold up in storage
η Efficiency [%] 40 Estimated
Fdemand Molar demand [mol s−1 ] 20 Estimated
25
Chapter 4
Integrated Scheduling and Control

Using Reduced Order Models
The framework includes several steps to establish the optimal schedule, see Figure 4.1. Since no historical
dataset was available for this process, the dynamic model established in Chapter 3 was used to simulate
the historical process data. First, the optimal steady state process operation was determined at several
production levels (from +20% of the nominal operation to −20% of the nominal operation). Then the
behavior of process variables during transitions between each production level were determined using
dynamic optimization. The reduced order model was then identified from the transition data. The
reduced order model was then used to find the optimal schedule of the production rate and liquid
inventory management. Finally, the optimal schedule was validated using the full dynamic model.
Training data
Steady state optimization Dynamic optimization
Reduced order model
Model identification
Scheduling
Dynamic optimization
Verification
Figure 4.1: Scheduling procedure
4.1 Optimal Production Levels

The nitrogen demand was assumed to be constant at 20 mol s−1 throughout the horizon. With the
storage tank, the production rate was free to deviate from the nominal value while still meeting the
demand. The optimal operating conditions at different production levels in the range of ±20% of the
nominal production level (20 mol s−1 ) was found using steady state optimization. The decision variables
were the rate of air feed to the plant, the fraction of the air stream liquefied in the multistream heat
exchanger, and the column reflux ratio. The objective was to minimize the plant’s energy usage at each
26
production level. Several constraints were imposed to ensure feasible operation at each production levels.
The solution to the optimization problem was found using gPROMS 4.0 [37] and the NLPSQP solver.
Decision Variables and Objective Function

The decision variables in the optimization were the rate of air fed to the process, the column reflux ratio,
and the fraction of air liquefied in the multistream heat exchanger. The optimal values of each decision
variable was obtained at every production level. The pressure exiting the compressor was held constant
at 5.7 bar and the pressure exiting the PHX turbine was set to 1.3 bar. The reboiler pressure was fixed
at 2.5 bar, and the nominal production rate was set to 20 mol s−1 equivalent of 50 Tons day−1 . The
optimal steady states were determined at several production levels which are listed in Table 4.1.
The optimal steady state operation seeks to minimize energy consumption. The objective function
imposed was:
φ(x)ss = (ω(Fdemand λ − Fproduction ))2 + (θWtotal )2 (4.1)
where ω and θ are penalty factors. λ is the deviation of production rate from the nominal value (Fdemand )
as given in Table 4.1, and Fproduction is the product flow rate. Wtotal is the total energy consumption in
the plant (not including the liquefier).
Table 4.1: Production levels used in the steady state optimization
Demand [mol s−1 ] λ

24 0.20
23 0.15
22 0.10
21 0.05
20 0
19 -0.05
18 -0.10
17 -0.15
16 -0.20
Heat Exchanger Constraints

As discussed in section 3.3.3, the heat exchanger was divided into two zones, where only sensible heat
was removed from the incoming air stream in zone 1, and the inlet air was condensed in zone 2. To
ensure that no phase transformations occurred in zone 1, the air feed pressure at the outlet of zone 1
was constrained to be in the vapor phase:
pdew − p1,4,50 ≥ 0 (4.2)
where pdew is the dew point pressure at the outlet of zone 1 [bar] and p1,4,50 is the pressure of the air
feed stream at the outlet of zone 1 [bar].
Additionally, a constraint was set to ensure that the outlet of the air feed stream from zone 2 was in
the liquid phase:
p2,4,1 − pbubble ≥ 0 (4.3)
where p2,4,1 is the pressure at the outlet of zone 2 [bar] and pbubble is the corresponding bubble point
pressure at the outlet of zone 2 [bar].
The dew point pressure was calculated as:
X
pdew = pdew,j xj (4.4)
j
where pdew is the dew point pressure of the mixture [bar], pdew,j is the dew point pressure of component
j [bar] calculated using equation (3.9), and xj is the liquid molar fraction of component j in the mixture.
27
Likewise, the bubble point pressure was calculated as:
X
pbubble = pbubble,j xj (4.5)
j
where pbubble is the bubble point pressure of the mixture [bar], pbubble,j is the bubble point pressure of
component j [bar] calculated using equation (3.9), and xj is the liquid molar fraction of component j in
the mixture.
The nitrogen product stream was at high pressure, and to fully utilize the available heat, the stream
was withdrawn from the heat exchanger and expanded in a turbine to produce electricity, see Figure 3.3.
To ensure that the nitrogen exiting the turbine was in the gas phase, the pressure was constrained to be
below the corresponding saturation pressure for nitrogen:
psat
N2 − pout ≥ 0 (4.6)
where psat
N2 is the saturation pressure for nitrogen [bar] calculated using equation (3.9) and pout is the
pressure exiting the turbine [bar].
Column Constraints
The flooding capacity constraint had to be met in the distillation column. At 100 % flooding, liquid is
excessively carried over to the adjacent stage due to a high vapor flow rate. This typically results in
drastically reduced stage efficiencies and high pressure drops. The flooding fraction is typically set at
80-90 % during nominal production [48].
The flooding constraint was adapted from the work by Coulson and Richardson [48] and Cao [38].
To ensure that no flooding occurred, the flooding fraction was constrained to be below 1:
ui
αif looding = <1 (4.7)
uf looding,i
where αif looding is the fractional flooding at stage i, ui is the vapor velocity at stage i [m s−1 ], and
uf looding,i is the flooding velocity at stage i [m s−1 ]. ui is calculated by:
Vi
ui = (4.8)
ρV AActive
where Vi is the molar vapor flow at stage i [mol s−1 ], ρV is the vapor molar density of the mixture at
stage i [mol m−3 ] calculated using equation (3.35), and AActive is the active area of stage i [m2 ].
The flooding velocity can be calculated as [38]:
s
ρL ML − ρV MV
uf looding,n = K1 (4.9)
ρV MV
where K1 is a capacity constant and ρL is the liquid density [mol m−3 ] obtained from equation (3.42).
ML is the molar weight of the liquid [kg mol−1 ] and MV is the molar weight of the vapor [kg mol−1 ].
K1 was calculated as: σn
K1 = K2 −
(4.10)
20 ∗ 10 3
where σn is the surface tension of the liquid mixture [N m−1 ] and K2 is a capacity factor calculated as:
K2 = (0.0744Ts + 0.0117) log F −1 + 0.0304Ts + 0.0153 (4.11)
where Ts is the tray spacing [m] and F is a flow parameter. The tray spacing was estimated from work
done by Thorogood et al. [49] and F was calculated as:
s
GL ρV MV
F = Vi (4.12)
Gi ρL ML
28
−1
where GL i is the mass flow rate of liquid at stage i [kg s ] and GVi is the mass flow rate of vapor at stage
−1
i [kg s ]. F was assigned to 0.1 if it was calculated to be lower than 0.1 [38].
The surface tension of the liquid mixture was estimated from the surface tension of the pure compo-
nents:
X 3
r
σm = xj σjr (4.13)
j
where σm is the surface tension of the mixture, σj is the surface tension of the pure component j
[mN m−1 ] and r is an adjustable parameter set to 0.25 [50]. The surface tension of the pure component
was calculated using an empirical relationship given by Kai et al. [50]:
Ti p
σj = σ0 (1 − ) (4.14)
Tc
where σ0 and p are parameters for each component j, Ti is the temperature at stage i, and Tc is the
critical temperature of component j [50].
Reboiler and Condenser Constraints

The product of the overall heat transfer coefficient and heat transfer area, U A, in the integrated re-
boiler/condenser changed with liquid level in the reboiler. A constraint was therefore enforce to ensure
that the heat flow was always from the condenser to the reboiler. The constraint was implemented by
ensuring that U A was positive:
UA ≥ 0 (4.15)
Furthermore, the molar holdup in the reboiler was constrained to not deplete over time. This was
enforced by ensuring that the liquid flow into the reboiler was greater than the vapor flow out of the
reboiler:
L30 − Vreboiler ≥ 0 (4.16)
where L30 is the molar flow of liquid into the reboiler [mol s−1 ] and Vreboiler is the molar flow of vapor
out from the reboiler [mol s−1 ]. The hold up was controlled by a level controller (equation (3.23)).
Finally, the maximum allowed impurity in the nitrogen product stream is 1000 ppm. To ensure that
this constraint was not violated during transitions, the steady state impurity is constrained below 500
ppm:
Impurity = (1 − yproduct )106 (4.17a)
Impurity ≤ 500ppm (4.17b)
where yprodcut is the mole fraction of nitrogen in the product stream.
4.2 Dynamic Optimization to Find Optimal Transitions

Dynamic optimization was used to simulate historical data of the plant transitioning between different
production levels. Optimal trajectories for the manipulated variables are determined which ensure a
rapid transition to the new steady state, and that no operational constraints are violated during the
transitions. This dataset was used to determine an approximate reduced order dynamic model of the
process. The solution to the optimization problem was found using gPROMS 4.0 [37] and the NLPSQP
solver.
Constraints
The constraints of the steady state optimization must hold during the transitions as well (equation (4.2)
- (4.7), (4.15) - (4.17)). The maximum allowable impurity during transitions was set to 1000 ppm. To
ensure that these constraints are satisfied throughout the time horizon, they were transformed into path
constraints. An example of this transformation for the dew point pressure constraint, equation (4.2), is
given by:
29
Z tf
max(0, p1,4,50 − pdew ) ≤ (4.18)
t0
where is a small positive parameter. These path constraints are implemented in gPROMS as endpoint
inequality constraints.
Discretization
To simulate the control moves, a piecewise linear control profile is used, see Figure A.1. To obtain the
piecewise linear control profile, the rate of change for each manipulated variable in each control interval
was determined:
dui,j
= ai,j (4.19)
dt
where ui corresponds to manipulated variable i (see Table 4.2) and ai,j is the rate of change of variable i
in control interval j. Thus, the optimization decision variables were actually piecewise constants, which
correspond to the slope of the manipulated variables in each control interval.
Table 4.2: Manipulated variables
Manipulated Variable Explanation

F1,4 Feed Flow
rgasdraw Fraction of product entering the condenser
kwithdrawn Fraction liquid air entering the column
five control intervals were established for each transition with an extra interval at the end to ensure
that steady state is reached (fix ai,j = 0 in the last interval). An additional constraint was imposed to
ensure that the optimal steady state values of the manipulated variables were reached at the end of the
time horizon [38]:
ui 2
(1 − ) ≤ (4.20)
uop
i,j
where uop
i are the optimal values of manipulated variable i at production rate j and is a very small
parameter.
Objective function
The objective of the dynamic optimization was to minimize the transition time between production
levels. The objective function was formulated as [38]:
" #2
Z tf
Fproduction (t)
φDO = 1− ∗ dt + t2f (4.21)
t0 Fproduction
∗
where Fproduction is the production rate at the end of the transition.
4.3 Reduced Order Model Identification

An accurate reduced order model used in an optimal scheduling formulation should preserve most of
the dynamic behavior of relevant process variables. It should also contain fewer state variables than a
full dynamic model. Here, the reduced order model was identified from the optimal transitions dataset
but any kind of operating data could be used. For this system, transfer function models were able to
capture the dynamic behavior accurately for most process outputs, and a nonlinear model was necessary
to capture the dynamic behavior of one output.
30
Relevant Process Variables
The desired production rate was the scheduling decision, and was therefore the input to the reduced
order model. The outputs were the process variables relevant to scheduling decisions. These include the
actual production rate and energy use, as well as the constrained operating variables. The constrained
variables must be included to ensure that the production rate changes are physically realizable. The
energy use was directly proportional to production rate, so only one model must be identified for these
two variables. The constrained operating variables included the product impurity, column flooding, the
phase of the inlet stream exiting zone 1 of the PHX, the phase of the inlet stream exiting zone 2 of the
PHX, the product of the overall heat transfer coefficient and area in the reboiler, and the reboiler liquid
level. The output variables and their corresponding identified model order are presented in Table 4.3.
Table 4.3: Outputs in the scale bridging model and model order
Output Model Order

Impurity 7
Flooding 8
Dew pressure in PHX Non-Linear
Bubble pressure in PHX 10
Production rate 8
UA in reboiler 8
Reboiler liquid level 9
Transfer Function model

The optimal transition dataset was used to identify the reduced order models. All but one important
output variable were accurately approximated by transfer functions. The model input was the desired
production flow rate (which will be determined by the schedule) and the outputs were the relevant process
variables discussed in the previous section (presented in Table 4.3). It was observed that the manipulated
variables, u, tended to follow similar optimal trajectories during the production level transitions. Thus,
to obtain more information about the dynamics of the process and to better fit the output variables,
the reduced system was modeled as a network of transfer functions which include the dynamic behavior
of the manipulated variables during the transitions. The outputs (y) were functions of the manipulated
variables (u) which were functions of the desired production rate set point (Fsp ) (depicted in Figure
4.2). First, transfer functions relating the production level setpoint and the manipulated variables were
identified:
U = Gu Fsp (4.22)
where Gu are transfer functions relating the manipulated variables to the production rate set point.
Next, transfer functions were identified which relate the manipulated variables to the relevant output
variables:
Y = Gy U (4.23)
where Gy are transfer functions which relate the 3 manipulated variables to each output. The overall
transfer function was then given by:
G = Gu Gy (4.24)
where G directly relates the production rate setpoint, Fsp , to the outputs. The complete network of
1x1 transfer functions is given in Figure 4.2 where Gi corresponds to a second order single input single
output transfer function. The overall transfer function was given more specifically by:
Gj = G1,j + G2,j G5,j + G3,j G6,j + G4,j G7,j (4.25)
Due to the fact that each individual transfer function for each output j, Gi,j , was second order, the
overall transfer function Gj was 14th order for each output j.
The parameters for the transfer functions were found using the System Identification Tool box in
MATLAB [51]. It was determined that the 14th order model had poor conditioning (i.e., the ratio of
31
G1,j
+
+
u1
ysp G2,j G5,j + y
+
u2
G3,j G6,j
u3
G4,j G7,j
Figure 4.2: Block diagram of the transfer function.
the largest and smallest eigenvalue is large). This causes difficulties for the simulator and results in very
long simulation times [4]. Therefore, the model is further reduced to eliminate the nonessential poles.
The balanced Gramian state space realization are obtained using the MATLAB [51] function balreal [52].
The diagonal entries of the joint Gramian, g, are evaluated and poles corresponding to g < 0.002 were
removed. The resulting model was of varying order between seven and ten, see Table 4.3.
Non-Linear Hammerstein Wiener model

The model for the dew pressure ratio exiting zone 1 of the PHX was not adequately described by a linear
transfer function. A Hammerstein-Wiener model was used to capture the nonlinear behavior. A second
order polynomial transforms both the inputs and outputs of the transfer function (third order) model:
2
u = AFsp + BFsp + C (4.26a)
0
y = Gu (4.26b)
y = ay 02 + by 0 + c (4.26c)
where G, the third order transfer function, and the parameters A, B, C, a, b, c were found using System
Identification Tool box in MATLAB [51].
The identified reduced models were transformed into state space models. The complete reduced order
process model also included a mass balance over the storage tank, equation (3.58), and an equation
ensuring that the demand is satisfied, equation (3.57). Furthermore, the energy requirement given by:
E(t)total = W (t)compressor + W (t)turbine,1 + W (t)turbine,2 + W (t)liquef ier (4.27)

Was calculated assuming that the compressor and turbine power requirement are directly proportional to
the inlet and production flowrates. W (t)compressor is the compressor power demand calculated by equa-
tion (3.43), W (t)turbine is the turbine power production calculated by equation (3.50), and W (t)liquef ier
is the power demand of the liquefier. Hpoly in equation (3.43) and equation (3.50) was assumed to be
constant. W (t)liquef ier was related to the liquefaction rate by:
W (t)liquef ier = F HL (4.28)

−1
where F is the molar flow being liquefied and HL is the polytropic head of the liquefier [J mol ], which
is assumed to be constant.
32
4.4 Scheduling Using the Reduced Order Model
The reduced order model was used to obtain the optimal schedule which minimizes the operating cost
of the plant while meeting product demand throughout the horizon. The scheduling problem setup was
formulated as a dynamic optimization problem to determine the optimal production rate set point. The
model was implemented and optimized using gPROMS 4.0.0 [37] and the NLPSQP solver.
Objective function and cost schedule

The objective function in this case was to minimize the operating cost throughout the scheduling horizon:
Z tf
φSO = E(t)total price(t)dt (4.29)
0
where Etotal is the power usage of the plant at time t [MW h] and price is the hourly electricity price
[EUR/MWh].
The electricity price schedule was taken from the Nord Pool Spot market for UK [53]. N2EX Day
Ahead Auction prices from April 11 - 13 2015 are presented in Figure 4.3. According to Ierapetritou et
al. [11] the price can be known with reasonable certainty over a 3 day time horizon. Predicted prices on
a longer time horizon are more uncertain and require optimal scheduling under uncertainty [54].
110
100
90
Electricity Price [Eur/MWh]
80
70
60
50
40
30
0 4 8 12 16 20 24 28 32 36 40 44 48 52 56 60 64 68 72
Time [hour]
Figure 4.3: Price schedule for April 11 - 13 2015
Decision variables and constraints

The constraints in the scheduling problem were equivalent to those in the dynamic optimization. These
included: i) the feed stream pressure exiting zone 1 of the PHX must be below the dew point pressure,
ii) the pressure exiting zone 2 of the PHX must be above the bubble point pressure, iii) the impurity
must never exceed 1000 ppm, iv) the column should never reach 100% flooding, v) the heat transfer
in the reboiler/condenser must be in the right direction, and vi) the reboiler liquid must not disappear
(equation (4.2), (4.3), (4.17), (4.7), (4.15) and (4.16), respectively were transformed to path constraints).
Moreover, to ensure physical realizability the impurity was not allowed to decrease below 0:
33
Impurity ≥ 0 (4.30)
Constraints are also imposed on the storage tank inventory. At the initial state the storage tank was
assumed to be 50 % full and at the end of the time horizon the storage tank was constrained to be at
least the same.
M (t0 )storage
∗ = 0.5 (4.31)
Mstorage
M (tf )storage
∗ ≥ 0.5 (4.32)
Mstorage
where M (t0 )storage is the liquid level in the storage tank at the beginning of the time horizon, M (tf )storage
∗
is the liquid level in the storage tank at the end of the time horizon and Mstorage is the maximum storage
tank capacity. Furthermore, constraints were implemented to ensure that the tank was not emptied or
filled beyond its capacity:
M (t)storage
0≤ ∗ ≤1 (4.33)
Mstorage
where M (t)storage is the liquid level in the storage tank at time t.
The scheduling decision variable, the production rate, was allowed to vary every hour by up to ±20%
of the nominal production rate. Furthermore, constraints were imposed on the rate of change of the
decision variable. In one hour, the production rate set point could not change by more then 3 mol s−1 :
dFproduction
≤3 (4.34)
dt
4.5 Implementation of Schedule on the Full Dynamic Model

The optimal schedule found with the reduced order model was implemented on the full dynamic model,
established in section 3.3, in order to evaluate the effectiveness and accuracy of the proposed scheduling
method.
Similar to the problem setup described in section 4.2, the optimal transitions between the determined
production rates were determined using dynamic optimization. The constraints were equivalent to those
presented earlier, equation (4.2) - (4.7), (4.15) - (4.17), but the objective function was slightly different
[38]:
" #2
Z tf
Fproduction (t)
min φIO = 1− ∗ dt (4.35)
a,δ t0 Fproduction
The decision variables were equivalent to the dynamic optimization, equation (4.19). The time horizon
was divided into 6 (variable length) control intervals for every hour. First, 5 control intervals were used
where ai,j was free to be optimized, and the last where ai,j = 0. At the sixth interval in every hour,
constraint (4.20) was imposed to ensure that the optimal values of ui,j were reached. To find the optimal
control values for each level, a polynomial was fit to interpolate the optimal values of ui,j found at each
discrete operating level in the steady state optimization, see Table 4.1. A first order polynomial was
used for the feed flow and a second order polynomial fit rgasdraw and kwithdraw .
34
Chapter 5
Results and Discussion
The results from the dynamic optimization, model identification and schedule optimization are presented
in this chapter, along with a simulation of the full dynamic model when subjected to the optimal schedule.
The results are normalized except when otherwise noted.
5.1 Model identification

The purpose of the dynamic optimization is to generate a set of simulated transition data for the ASU.
No real process data was available for this work but the framework presented can easily be adapted to
utilize true plant transition data. A set of production setpoint transitions was selected to create the
training dataset that traverses the range of possible production rates, as can be seen in Figure 5.1. The
dynamic optimizations found the optimal manipulated variable trajectories for each transition (Table
4.2). The relevant output variable trajectories (Table 4.3) are presented in Figure 5.2 except for the PHX
zone 1 exit pressure which is presented in Figure 5.3. The Normalized Mean Square Error (NRMSE) is
given between the identified model and the data and it is defined as [51]:

||xref − x||
N RM SE = 100 1 − (5.1)
||xref − mean(xref )||
where xref is the reference data vector and x is the test data vector.
35
24
Production rate set point [mol/s]

23
22
21
20
19
18
17
16
0 10 20 30 40 50 60 70 80
Time [hour]
Figure 5.1: Step changes in production set point considered during the dynamic optimization
As seen in Figure 5.2a, 5.2b, and 5.3 the impurity, the column flooding and the ratio between the dew
pressure and the pressure exiting zone 1 of the PHX (denoted dew pressure ratio henceforth) are all close
to or temporarily violate their corresponding constraints during the transition horizon. The impurity
level, which has a very slow dynamic behavior, is the most challenging constraint to satisfy during the
transitions. Moreover, the ratio between the pressure of the feed stream exiting zone 2 of the PHX and
the bubble pressure, see Figure 5.2c, (henceforth denoted bubble pressure ratio) also approaches the
constraint bounds. The heat transfer constraint, and the reboiler liquid level are not close to the bounds
throughout the trajectory, see Figure 5.2e and 5.2f.
The dashed line in Figure 5.2 shows the outputs of the reduced model, equation (4.24). The NRMSE
between the training data outputs and the reduced model prediction is also given. For production,
flooding and heat transfer the identified reduced model predictions fit the data very well with NRMSE
values between 75 % - 88 %. Especially the model for the production rate, figure 5.2d, has a good fit
with an NRMSE value of 88 %. The NRMSE between the bubble pressure ratio training data and the
corresponding reduced order model prediction was relatively low. At 25 hours and 60 hours, a significant
offset is observed between the data and the model prediction. However, most of the transitions are
accurately predicted, and the transfer function model for the bubble pressure ratio captures the dynamic
behavior of the training data acceptably. The reduced transfer function model for the impurity level
did not have a good fit (NRMSE of 34 %), due to a slight mismatch in the gain, but the time constant
appears to mimic the dynamic behavior of the process. It was determined that the fit was acceptable,
and if the schedule suggests a transition which may lead to a breach of the impurity constraint, the
underlying control system should be able to prevent these occurrences.
The reduced order transfer function model of the dew pressure ratio did not have a sufficient fit with
a NRMSE value of -74 %. In Figure 5.3, the model predicts constraint violations at several points when
the training data does not violate the constraints. A nonlinear Hammerstein-Wiener model was therefore
identified to better predict the dew point pressure. As can be seen in Figure 5.3 the nonlinear reduced
order model for the dew pressure ratio gave a significantly better prediction than the transfer function
model.
36
1
2
Data 0.8 Data
Transfer Function model Transfer Function model
1.5 Constraint level 0.6 Constraint Level
NRMSE: 34.1463 % NRMSE: 75.9094 %
0.4
1 0.2
Flooding
Impurity
0
0.5
−0.2
0 −0.4
−0.6
−0.5
−0.8
−1 −1
0 10 20 30 40 50 60 70 80 0 10 20 30 40 50 60 70 80
Time [hour] Time [hour]
(a) Impurity (b) Flooding

1.5
4
1
2
0.5
Bubble pressure
Production
0
0
-2
−0.5
-4 Data
Constraint Level −1
NRMSE: 88.1496 %
NRMSE: 38.7112 % Data
-6
−1.5
0 10 20 30 40 50 60 70 80 0 10 20 30 40 50 60 70 80
(c) Bubbel pressure (d) Production

4 2.5
Data
2 Transfer Function model
3 NRMSE: −59.285 %
1.5
2 1
Reboiler Liquid level
0.5
1
0
UA
−0.5
0
−1
−1 −1.5
−2
−2 NRMSE: 83.2343 %
Data −2.5
−3 −3
0 10 20 30 40 50 60 70 80 0 10 20 30 40 50 60 70 80
(e) UA (f ) Reboiler liquid level
Figure 5.2: Data and transfer function models for step change responses.
37
4
Data
Transfer Function model, NRMSE: −74.6451 %
3.5 Non−linear model, NRMSE: 61.3946 %
Constraint level
3
2.5
Dew pressure
1.5
0.5
−0.5
−1
0 10 20 30 40 50 60 70 80
Time [hour]
Figure 5.3: Dew pressure data compared to transfer function model and non-linear model
Overall, the transition dataset obtained through dynamic optimization resulted in highly nonlinear
behavior of the model outputs. This is not surprising given the complex, highly coupled, and nonlinear
behavior of the equations that describe an ASU. In light of this fact, it is very encouraging that the
identified reduced order models obtained using transfer functions and Hammerstein-Wiener models were
able to predict the important dynamic behaviors of the process variables during open-loop production
rate transitions. The fits are not perfect, but the models are suitable for optimal scheduling calculations
with a significantly smaller model. Additionally, a control system would be utilized to ensure that the
constraints are satisfied during the optimal schedule.
5.2 Scheduling Using the Reduced Order Model

Figure 5.4 shows the optimal schedule found using the reduced order models and the optimization problem
formulated in section 4.4 over a 72 hour horizon (1 corresponds to a 20 % increase from the nominal
production rate of 20 mol s−1 and -1 corresponds to a 20 % decrease in production rate).
The total operating cost for 72 hours with a fixed production rate at the nominal value is e31,266.
With the optimal production rate schedule the cost is e30,466, which corresponds to a 2.6 % cost
reduction over the 3 days considered. In accordance with the price schedule (Figure 4.3) the production
rate is ramped up when prices are low and ramped down when prices are high. Furthermore, the storage
tank was periodically filled and emptied during the 72 hours, as can be seen in Figure 5.5, where 1
corresponds to a filled tank and 0 corresponds to an empty tank. This shows that the plant was fully
utilizing the storage tank capacity in order to minimize cost, as it was emptied during high price periods,
and filled during low price periods. It also shows that the operating cost savings resulting from a variable
production rate operation depends on the size of the storage tank. The smaller the storage tank, the
more rapidly the tank is filled and emptied, and the less time that the plant can produce above or below
the product demand rate during low and high price scenarios, respectively.
38
1
0.8
Production set point 0.6
0.4
0.2
−0.2
−0.4
−0.6
−0.8
−1
0 10 20 30 40 50 60 70 80
Time [hour]
Figure 5.4: Schedule for production rate set point
0.9
0.8
Hold Up Storage Tank
0.7
0.6
0.5
0.4
0.3
0.2
0.1
0
0 10 20 30 40 50 60 70 80
Time [hour]
Figure 5.5: Hold up in storage unit
The relevant process variables (reduced order model outputs) are presented in Figure 5.6. Both
the flooding and bubble pressure ratio reduced model predictions do not approach their constraints,
39
see Figure 5.6c and Figure 5.6d. The reduced model prediction of the dew pressure ratio was close to
the constraint violation throughout the schedule horizon, see Figure 5.6b. The impurity reduced model
prediction of the impurity level, Figure 5.6a, reaches its constraint bound on several occasions. Due to
the fact that the impurity reduced model prediction was the least accurate, it is important that a robust
controller is implemented to ensure that the actual system does not violate the constraint throughout
the horizon. The other reduced model output predictions are presented in the appendix.
1.5
Impurity 3 Dew Pressure
Constraint level
Constraint level
1 2.5
Dew Presure
0.5
Impurity
1.5
0 1
0.5
−0.5
0
−1 −0.5
0 20 40 60 80 0 20 40 60 80
(a) Impurity (b) Dew pressure

1 3
Flooding
Constraint level 2
1
0.5
0
Bubbel Pressure
−1
Flooding
0 −2
−3 Bubble pressure
−4 Constraint level
−0.5
−5
−6
−1 −7
0 20 40 60 80 0 20 40 60 80
(c) Flooding (d) Bubble pressure
Figure 5.6: Expected constraint and variable profiles from the dynamic scheduling using the reduced order model
for 72 hours
5.3 Implementation of schedule in the full dynamic model

Figure 5.7 shows the production rate determined by the full dynamic model and the prediction given by
the reduced order model when subject to the optimal production rate schedule. The NRMSE between the
reduced order model prediction and the full dynamic model is also given. The production rate set point
changes correspond to the first 24 hours of the optimal schedule found using the reduced order model.
Only the first 24 hours were able to be implemented because the dynamic optimization over a longer
time horizon requires several weeks of computation time to arrive at a solution. As seen in the figure,
the production rate prediction closely tracked the rate determined by the full dynamic model. The slight
discrepancies occurred mainly at low production rates which causes divergence of the storage tank hold
up, as can be seen in Figure 5.8. Recall that the training dataset and the testing data were obtained using
separate dynamic optimizations (i.e., the manipulated variable trajectories are not identical during the
transitions). Nevertheless, the electricity cost calculated with the full dynamic model and that predicted
by the reduced order model was about the same: e11,032 for the dynamic model and e11,039 for the
reduced order model over the first day of the schedule. The storage tank capacity is important to take
into account. As can be seen in Figure 5.8 the full dynamic model shows that the storage tank will
empty beyond its capacity. A buffer should be implemented in the optimal scheduling calculations (e.g.,
the tank cannot empty beyond 10% of the capacity throughout the horizon).
40
0.8 Production set point
Full dynamic model
0.6 Reduced model
NRMSE: 82.753 %
0.4
0.2
Production
−0.2
−0.4
−0.6
−0.8
−1
0 2 4 6 8 10 12 14 16 18 20 22 24 26
Time [hour]
Figure 5.7: Production of the plant using the dynamic and reduced order model for the first 24 hours of the
schedule found using the reduced order model. NRMSE between the reduced order model and the full dynamic
model is also given
0.9
0.8
0.7
Hold Up Storage Tank
0.6 NRMSE: 91.6234 %

Full dynamic model
0.5 Reduced model
Constraint level
0.4
0.3
0.2
0.1
−0.1
0 2 4 6 8 10 12 14 16 18 20 22 24
Time [hour]
Figure 5.8: Hold up in storage during 24 hours using the dynamic and reduced model
41
The impurity level is shown in Figure 5.9a for the full dynamic model and the reduced order model
prediction. While the NRMSE between the reduced model prediction and the full model is 36 %,
the model accurately predicts the impurity spikes, and captures similar dynamic features of the full
dynamic model. The reduced model prediction of the impurity during the optimal 3 day schedule
(Figure 5.6a) suggests that the impurity level would reach its bounds on several occasions. The impurity
level determined using the full dynamic model over the first 24 hours indicates that these bounds are in
fact reached as predicted by the reduced order model.
In Figure 5.9b, the dew pressure ratio determined by the full dynamic model and the reduced model
prediction is shown over the first 24 hours of the optimal schedule. The reduced order model accurately
predicts the time periods in which the dew pressure ratio is close to the constraint bounds. Other
model outputs including the flooding fraction, the bubble pressure ratio, and the heat transfer constraint
display a tight fit between the reduced order model predictions and the full dynamic model outputs (see
Figures 5.9c, 5.9d, and 5.9e). The reboiler liquid level reduced model is unable to satisfactorily predict
the dynamic behavior determined by the full model, however, the variable is not close to the constraint
bound throughout the horizon (see Figure 5.9f).
The results indicate that the identified reduced order models are capable of finding an optimal feasible
production rate schedule that accurately predicts the transient behavior of relevant process variables. The
reduced order model consisted of 525 equations and 67 differential variables compared to the full dynamic
model which consisted of 5079 equations and 437 differential variables. The resulting optimization
problem formulation is tractable and can be solved in a relatively short time frame (14861 CPU seconds)
compared to an equivalent optimization formulation with a full dynamic model (which would probably
take weeks to solve with the most advanced optimization solvers). The results should be compared to
those determined using a static scheduling approach, where the production rate setpoint schedule is
determined under the assumption that some transitions are infeasible. While this approach is much
simpler, it is unlikely that the constrained variables will pushed to their bounds, or that the optimal
schedule is actually feasible. Consequently, it is expected that the optimal schedule determined using a
dynamic model will find a better result.
42
1 Full dynamic model
3.5
Reduced model
0.8 Constraint level
3 NRMSE: 39.4096 %
0.6
2.5
0.4
Dew Pressure
0.2 2
Impurity
0
1.5
−0.2
1
−0.4
0.5
−0.6
NRMSE: 36.0379 %
−0.8 Full dynamic model 0
Reduced model
Constraint level
−1
−0.5
0 2 4 6 8 10 12 14 16 18 20 22 24 0 2 4 6 8 10 12 14 16 18 20 22 24
(a) Impurity (b) Dew pressure

1
3
0.8
2
0.6
1
0.4
0
Bubbel Pressure
0.2
Flooding
−1
0
−2
−0.2
−3
−0.4
−4
−0.6
NRMSE: 66.3762 % −5 NRMSE: 36.5008 %
Full dynamic model Full dynamic model
−0.8 Reduced model −6 Reduced model
Constraint level Constraint level
−1
−7
0 2 4 6 8 10 12 14 16 18 20 22 24 0 2 4 6 8 10 12 14 16 18 20 22 24
(c) Flooding (d) Bubble pressure

2
2 NRMSE: 58.4499 %
Full dynamic model 1.5
Reduced model
1.5
1
1 0.5
Reboiler Liquid Level
0.5 0
−0.5
UA
0
−1
−0.5
−1.5
−1
−2
−1.5 −2.5 NRMSE: −34.2085 %
Full dynamic model
−2 −3 Reduced model
0 2 4 6 8 10 12 14 16 18 20 22 24 0 2 4 6 8 10 12 14 16 18 20 22 24
(e) UA (f ) Reboiler liquid level
Figure 5.9: Variables and constraints using the full dynamic model and the reduced order model.
43
Chapter 6
Conclusions and Future Work
A novel framework for making optimal scheduling decisions while incorporating process dynamics is
described in this thesis. The concept is to establish a reduced order model that predicts the dynamic
behavior of relevant process variables, e.g., those involved in operational constraints or process outputs,
when subject to changes in the production setpoint schedule. The optimal schedule is then determined
using the reduced order dynamic process model. This significantly differs from the conventional methods
using static models of full dynamic process models. The framework provides the benefit of using the
dynamic process knowledge to place process variables at their limits (if needed) in order to more rapidly
transition between production levels, while maintaining the process within its operational constraints.
The benefits of the framework are illustrated with an ASU that produces nitrogen at 99.9 % purity and
subjected to time-varying electricity prices. A detailed dynamic model of the ASU was developed which
included a cryogenic distillation column with reboiler/condenser, a PHX, a compressor, two turbines
and a liquefier, and consisted of 5079 equations with 437 differential variables. A simulated production
rate transition dataset was obtained using dynamic optimization on the dynamic model to determine the
optimal manipulated variable trajectories throughout the selected transitions. These data were used to
identify an accurate reduced order model consisting of 525 equations and 67 differential variables which
captures the dynamic behavior of important process variables (i.e., process outputs and constrained
variables). The reduced model utilizes transfer function models for all but one output and a nonlinear
Hammerstein-Wiener model for the the other output. An optimal production rate setpoint schedule
was determined for the ASU when subject to time-varying electricity prices with the use of a storage
tank. The optimization sought to minimize electricity cost throughout the time horizon by adjusting the
production rate setpoint while ensuring that the operational constraints were satisfied. As expected, the
optimal schedule utilized the stored liquid nitrogen to meet demand when electricity prices were high,
and filled the storage tank with liquid nitrogen when electricity prices were low. The optimal schedule
realized a 2.6 % savings compared to a flat production rate over the 3 day period considered. Finally,
the optimal schedule was implemented in the full dynamic model and again the optimal manipulated
variable transition trajectories were determined using dynamic optimization. The result showed that the
optimal schedule was feasible and that the predicted outputs from the reduced order model showed a
close match to their corresponding variables in the full model.
The limitation of the proposed framework is the ability to develop reduced order models that can
accurately predict the important process variables throughout production transitions. Given the complex,
highly coupled, and nonlinear nature of the equations that describe an ASU, it is encouraging that
accurate reduced order models can be obtained for the important process variables using relatively few
state variables. It was observed from the optimal trajectories that the manipulated variables followed
certain paths. When constructing the reduced order model it was advantageous to first related the input
variable to the manipulated variables with a transfer functions. The manipulated variables was thereafter
related to the output variables to obtain the overall transfer function. This structured approach gave
good results when identifying the reduced order model.
The framework presented in this thesis is generally applicable to optimal scheduling of processes
with variable product qualities and production rates. It is an alternative to the conventional scheduling
approaches which either use a static model with the process dynamics lumped into transition times, or a
full dynamic process model. In this work a full dynamic model of the plant was necessary to simulate the
44
transition data, but the framework can be applied when only historical transition data from the plant is
available (i.e., a full dynamic process model is not required if rich historical data is present). Finally, in
the presence of unexpected disturbances, the framework provides a quick and robust way to re-schedule
the plant operation on-line.
6.1 Future Work

Some suggestions for future work on this topic is presented in this section.
Comparing static scheduling to the reduced order model

A comparison with the conventional static scheduling approach should be performed in order to evaluate
the difference in the obtained results. In the static scheduling approach a table of allowable transitions
and other dynamic information is used to constrain the setpoint changes. Both the electricity cost
savings and the potential constraint violations should then be compared between the static schedule and
the reduced order model schedule.
Implementing a closed-loop control system

In this work, an open loop control system was used to obtain the transition dataset by using dynamic
optimization to determine the optimal manipulated variable trajectories. The resulting transition data
was highly nonlinear and non-smooth. In reality, a closed-loop control system would be in place to ensure
that some of the variables stay near their setpoints. The result could potentially be a more smooth and
linear dynamic behavior of the constrained process variables which are much more easily predicted with
a low order model. However, modeling challenges could arise if input saturation, non-smoothness, or
non-linear behavior exists in the data.
45
Bibliography
[1] M.M. El-Halwagi. Introduction to process integration. In Process Integration, volume 7 of Process
Systems Engineering, pages 1 – 20. Academic Press, 2006.
[2] M. Baldea and I. Harjunkoski. Integrated production scheduling and process control: A systematic
review. Computers & Chemical Engineering, 71:377–390, 2014.
[3] D.E. Shobrys and D.C. White. Planning, scheduling and control systems: why cannot they work
together. Computers & Chemical Engineering, 26(2):149 – 160, 2002.
[4] D.E. Seborg, D.A. Mellichamp, T.F. Edgar, and F.J. Doyle III. Process dynamics and control. John
Wiley & Sons, 2010.
[5] T.F. Edgar. Control and operations: when does controllability equal profitability? Computers &
Chemical Engineering, 29(1):41 – 49, 2004. {PSE} 2003.
[6] J. Du, J. Park, I. Harjunkoski, and M. Baldea. A time scale-bridging approach for integrating
production scheduling and process control. Computers & Chemical Engineering, 79(0):59 – 69,
2015.
[7] I. Harjunkoski, R. Nyström, and A. Horch. Integration of scheduling and control—theory or practice?
Computers & Chemical Engineering, 33(12):1909 – 1918, 2009. {FOCAPO} 2008 – Selected Papers
from the Fifth International Conference on Foundations of Computer-Aided Process Operations.
[8] J. Park, J. Du, I. Harjunkoski, and M. Baldea. Integration of scheduling and control using internal
coupling models. In J.J. Klemes, P.S. Varbanov, and P.Y. Liew, editors, 24TH European Symposium
on Computer Aided Process Engineering, PTS A and B, volume 33 of Computer-Aided Chemical
Engineering, pages 529–534, 2014.
[9] Y. Cao, C.L.E. Swartz, and M. Baldea. Design for dynamic performance: Application to an air
separation unit. In American Control Conference (ACC), 2011, pages 2683–2688, June 2011.
[10] M.H. Karwan and M.F. Keblis. Operations planning with real time pricing of a primary input.
Computers & operations research, 34(3):848–867, 2007.
[11] M.G. Ierapetritou, D. Wu, J. Vin, P. Sweeney, and M. Chigirinskiy. Cost minimization in an energy-
intensive plant using mathematical programming approaches. Industrial & engineering chemistry
research, 41(21):5262–5277, 2002.
[12] D.R. Vinson. Air separation control technology. Computers & Chemical Engineering, 30(10):1436–
1446, 2006.
[13] A.R. Smith and J. Klosek. A review of air separation technologies and their integration with energy
conversion processes. Fuel Processing Technology, 70(2):115 – 134, 2001.
[14] R. Huang, V.M. Zavala, and L.T. Biegler. Advanced step nonlinear model predictive control for air
separation units. Journal of Process Control, 19(4):678 – 685, 2009.
[15] T. Guo, J. Lu, W. Xiang, Y. Guo, and K. Wu. Dynamic modeling and control of the air separation
unit in an igcc power plant. In Sustainable Power Generation and Supply, 2009. SUPERGEN’09.
International Conference on, pages 1–7. IEEE, 2009.
46
[16] B. Roffel, B.H.L. Betlem, and J.A.F. De Ruijter. First principles dynamic modeling and multivari-
able control of a cryogenic distillation process. Computers & Chemical Engineering, 24(1):111–123,
2000.
[17] Y. Zhu, S. Legg, and C.D. Laird. Optimal design of cryogenic air separation columns under uncer-
tainty. Computers & Chemical Engineering, 34(9):1377 – 1384, 2010. Selected papers from the 7th
International Conference on the Foundations of Computer-Aided Process Design (FOCAPD, 2009,
Breckenridge, Colorado, USA.
[18] Energy Information Association. Manufacturing energy consumption survey: Total consumption of
electricity, 2010. http://www.eia.gov/consumption/manufacturing/data/2010/pdf/Table11_
1.pdf.
[19] R. Agrawal and T.F. Yee. Heat pumps for thermally linked distillation columns: An exercise for
argon production from air. Industrial & engineering chemistry research, 33(11):2717–2730, 1994.
[20] R. Agrawal, D.W. Woodward, K.A. Ludwig, and D.L. Bennett. Impact of low pressure drop struc-
tured packing on air distillation. In Institution of Chemical Engineers Symposium Series, volume
128, pages A125–A125. Hemsphere Publishing Corporation, 1992.
[21] M. Aneke and M. Wang. Potential for improving the energy efficiency of cryogenic air separation
unit (ASU) using binary heat recovery cycles . Applied Thermal Engineering, 81(0):223 – 231, 2015.
[22] J. Rizk, M. Nemer, and D. Clodic. A real column design exergy optimization of a cryogenic air
separation unit. Energy, 37(1):417 – 429, 2012. 7th Biennial International Workshop “Advances in
Energy Studies”.
[23] Yasuki Kansha, Akira Kishimoto, Tsuguhiko Nakagawa, and Atsushi Tsutsumi. A novel cryogenic
air separation process based on self-heat recuperation. Separation and Purification Technology,
77(3):389 – 396, 2011.
[24] W.F. Castle. Air separation and liquefaction: recent developments and prospects for the beginning
of the new millennium. International Journal of Refrigeration, 25(1):158 – 172, 2002.
[25] B. Keenan and K. Reuter. The design, experience and economic advantage of cryogenic turbo-
expanders with active magnetic bearings. In 3rd International Conference Kryogenica’94, 1994.
[26] R. Boehme, J.A.R. Parise, and R. Pitanga Marques. Simulation of multistream plate–fin heat
exchangers of an air separation unit. Cryogenics, 43(6):325 – 334, 2003.
[27] A.R. Smith, J. Klosek, and D.W. Woodward. Next-generation integration concepts for air separation
units and gas turbines. Journal of Engineering for Gas Turbines and Power, 119(2):298–304, 1996.
[28] B. Daryanian, R.D. Tabors, and R.E. Bohn. Optimal demand-side response to electricity spot prices
for storage-type customers. IEEE Trans. Power Syst.;(United States), 4(3), 1989.
[29] J.G. Roos and I.E. Lane. Industrial power demand response analysis for one-part real-time pricing.
Power Systems, IEEE Transactions on, 13(1):159–164, 1998.
[30] Y. Zhu, S. Legg, and C.D. Laird. A multiperiod nonlinear programming approach for operation of
air separation plants with variable power pricing. AIChE Journal, 57(9):2421–2430, 2011.
[31] S. Mitra, I.E. Grossmann, J.M. Pinto, and N. Arora. Optimal production planning under time-
sensitive electricity prices for continuous power-intensive processes. Computers & Chemical Engi-
neering, 38(0):171 – 184, 2012.
[32] J. Miller, W.L. Luyben, and S. Blouin. Economic incentive for intermittent operation of air sepa-
ration plants with variable power costs. Industrial & Engineering Chemistry Research, 47(4):1132–
1139, 2008.
[33] C. Zhongzhou, M.A. Henson, P. Belanger, and L. Megan. Nonlinear model predictive control of
high purity distillation columns for cryogenic air separation. Control Systems Technology, IEEE
Transactions on, 18(4):811–821, July 2010.
47
[34] G.Y. Zhu, M.A. Henson, and L. Megan. Low-order dynamic modeling of cryogenic distillation
columns based on nonlinear wave phenomenon. Separation and Purification Technology, 24(3):467
– 487, 2001.
[35] Y. Fu and X.G. Liu. Nonlinear control based on wave model of high-purity heat integrated air
separation column. Chemometrics and Intelligent Laboratory Systems, 146(0):1 – 9, 2015.
[36] R. Lopez-Negrete, F.J D’Amato, L.T Biegler, and A. Kumar. Fast nonlinear model predictive
control: Formulation and industrial process applications. Computers & Chemical Engineering,
51(0):55 – 64, 2013.
[37] Process Systems Enterprise. gPROMS 4.0.0. www.psenterprise.com/gproms, 1997 - 2014.
[38] Y. Cao. Design for dynamic performance: Application to an air separation unit, 2011.
[39] Aspen Technology, Inc. Aspen Plus. https://www.aspentech.com/products/aspen-plus.aspx,
1994 - 2015.
[40] C.L. Yaws. The Yaws Handbook of Vapor Pressure: Antoine coefficients. Elsevier Science, 2015.
[41] A. Harmens. Vapour-liquid equilibrium N2 -A2 -O2 for lower argon concentrations. Cryogenics,
10(5):406–409, 1970.
[42] W.M Haynes. CRC Handbook of Chemistry and Physics. CRC press, 2013.
[43] ALPEMA. The Standards of the Brazed Aluminum Plate-Fin Heat Exchangers Manufacturers’s
Association, 2000.
[44] D. Green and R. Perry. Perry’s Chemical Engineers’ Handbook, Eighth Edition. Chemical Engineers
Handbook. Mcgraw-hill, 2007.
[45] E.W. Lemmon and S.G. Penoncello. The surface tension of air and air component mixtures. In
P. Kittel, editor, Advances in Cryogenic Engineering, volume 39 of Advances in Cryogenic Engi-
neering, pages 1927–1934. Springer US, 1994.
[46] C.F. Spencer and S.B. Adler. A critical review of equations for predicting saturated liquid density.
Journal of Chemical & Engineering Data, 23(1):82–89, 1978.
[47] G. Venkatarathnam and K.D. Timmerhaus. Cryogenic mixed refrigerant processes. Springer, 2008.
[48] R.K. Sinnott and G. Towler. Chemical Engineering Design. Elsevier, 5th edition, 2009.
[49] R.M. Thorogood, D.L. Bennett, B.K. Dawson, A.L. Prentice, and R.J. Allam. Air separation process
using packed columns for oxygen and argon recovery, November 25 1992. Patent EP 0341854 B1.
[50] T. Kai, H. Nomoto, M. Deguchi, and T. Takahashi. Surface tension of ternary mixtures of nitrogen,
oxygen, and argon. Journal of Chemical & Engineering Data, 39(3):499–501, 1994.
[51] MathWorks, Inc. MATLAB Student R2014a. http://www.mathworks.com/products/matlab/,
1994 - 2014.
[52] A.J. Laub, M.T. Heath, C. Paige, and R. Ward. Computation of system balancing transforma-
tions and other applications of simultaneous diagonalization algorithms. Automatic Control, IEEE
Transactions on, 32(2):115–122, Feb 1987.
[53] Nord Pool Spot AS. N2EX Day Ahead Auction Prices. http://www.nordpoolspot.com/
Market-data1/N2EX/Auction-prices/UK/Hourly/?view=chart. [Online; acessed 15-May-2015].
[54] Q. Zhang, I.E. Grossmann, C.F. Heuberger, A. Sundaramoorthy, and J.M. Pinto. Air separation
with cryogenic energy storage: Optimal scheduling considering electric energy and reserve markets.
AIChE Journal, 61(5):1547–1558, 2015.
[55] S.S. Rao. Introduction to optimization. In Engineering Optimization, pages 1–62. John Wiley &
Sons, Inc., 2009.
48
[56] S. Skogestad. Self-optimizing control: the missing link between steady-state optimization and con-
trol. Computers & Chemical Engineering, 24(2):569–575, 2000.
[57] W.F. Feehery and P.I. Barton. Dynamic optimization with equality path constraints. Industrial &
Engineering Chemistry Research, 38(6):2350–2363, 1999.
[58] A. Cervantes and L.T. Biegler. Optimization strategies for dynamic systems optimization strategies
for dynamic systems. In C.A. Floudas and P.M. Pardalos, editors, Encyclopedia of Optimization,
pages 2847–2858. Springer US, 2009.
[59] V.S. Vassiliadis, R.W.H. Sargent, and C.C. Pantelides. Solution of a class of multistage dynamic
optimization problems. 1. problems without path constraints. Industrial & Engineering Chemistry
Research, 33(9):2111–2122, 1994.
[60] L.T. Biegler. Nonlinear programming strategies for dynamic chemical process optimization. Theo-
retical Foundations of Chemical Engineering, 48(5):541–554, 2014.
[61] S.S. Rao. Nonlinear programming iii: Constrained optimization techniques. In Engineering Opti-
mization, pages 380–491. John Wiley & Sons, Inc., 2009.
[62] J. Sjöberg, Q. Zhang, L. Ljung, A. Benveniste, B. Delyon, P.-Y. Glorennec, H. Hjalmarsson and, A.
Juditsky. Nonlinear black-box modeling in system identification: a unified overview. Automatica,
31(12):1691 – 1724, 1995.
[63] L. Ljung. System identification. In Wiley Encyclopedia of Electrical and Electronics Engineering.
John Wiley & Sons, Inc., 2001.
[64] D. Zill and M. Cullen. Differential Equations with Boundary-Value Problems. Cengage Learning,
2008.
[65] A.P. Wardle. Chapter 7 - Process Control . In J.F. Richardson D.G. Peacock, editor, Coulson &
Richardson’s Chemical Engineering Volume 3, pages 560 – 736. Butterworth-Heinemann, Oxford,
third edition edition, 1991.
[66] A. Wills, T.B. Schön, L. Ljung, and B. Ninness. Identification of Hammerstein–Wiener models .
Automatica, 49(1):70 – 81, 2013.
49
Appendix A
Concepts
A.1 Differential Alegbraic Equation System

The mathematical model of an ASU is typically presented in the form of a system of DAE that must
be solved numerically over a time horizon of interest. A simple example of a DAE system is presented
below:
x01 = 3x1 + x2 − 5y1 (A.1a)
x02 = x2 + 3y1 (A.1b)
y1 = x1 − 6x2 (A.1c)
where x are the differential variables and y are the algebraic variables. This system is classified as index
1 because the algebraic equation (A.1c), only needs to be differentiated once to transform the system
into a set of Ordinary Differential Equations (ODE):
y10 = x01 − 6x02 = 3x1 − 5x2 − 23y1 (A.2)
The system of ODE can be solved using several different numerical methods. Examples include the Euler
method, the trapezoidal method, the Taylor series method, the Runge-Kutta method, and the linear
multi-step method [38]. Additionally, several solvers (e.g., gPROMS [37]) are capable of integrating
index 1 DAE systems without transforming the algebraic equations. However, it is possible that higher
index DAE systems arise in process models. For example, if equation (A.1c) is replaced with:
0 = x1 − 6x2 (A.3)
The algebraic equation would need to be differentiated twice before the system can be transformed into
a set of ODE:
0 = x01 − 6x02 = 3x1 − 5x2 − 23y1 (A.4)
3x01 − 5x02
y10 = (A.5)
−23
This DAE system is of index 2. Such problems typically cannot be solved in this form. Hidden con-
straints in the algebraic equation, stability restrictions, and the difficulty in determining consistent initial
conditions are some of the problems with high index DAE. Dynamic models of distillation columns, for
example a cryogenic distillation column, are typically index 2 [14, 38]. To solve a high index DAE, the
problem needs to be reformulated into index 1 to reduce the order (see section 3.3.1 for an example of
index reduction). This can be costly depending on the size and complexity of the original DAE system.
Direct integration approaches without order reduction may be used, but the systems must have a certain
structure for this option to be viable [38].
50
A.2 Optimization
Broadly speaking, there are 2 classes of process optimization, steady state optimization and dynamic
optimization. Steady state optimization is used, for example, to find the optimal equipment design
and the optimal operating parameters, such as flow rates, temperatures, pressure, etc., to minimize the
capital and the operating cost of the nominal operation. Dynamic optimization is typically used to find
the optimal control moves between operating points over a time horizon to maximize the plant efficiency
and profit [4].
A.2.1 Optimization Problem Formulation

An optimization problem formulation includes an objective function, equality constraints, and inequality
constraints. The objective function is the figure of merit for the system, equality constraints typically
define the mathematical process model, and inequality constraints typically define operational or design
limitations. Variables in an optimization problem may be either continuous or discrete. Discrete variables
are often required to describe phase transitions or the presence (or lack of) equipment, and result in a
significantly more complex problem to solve. Solvers have been developed to handle different classes of
problems including linear and discrete (mixed integer) linear problems, as well as nonlinear and mixed
integer nonlinear problems [55].
A steady state optimization problem is formulated generally as [56]:
min φ(x, z, θ) (A.6a)

z
Subject to
F (x, z, θ) = 0 (A.6b)
g(x, z, θ) ≤ 0 (A.6c)
Likewise, the general form of a dynamic optimization problem is described as [57, 58]:
min φ(x(t), y(t), u(t), tf , θ) (A.7a)

u(t),tf
Subject to:
f (ẋ(t), x(t), y(t), u(t), t, θ) = 0 (A.7b)
h(x(t), y(t), u(t), t, θ) = 0 (A.7c)
g(x(t), y(t), u(t), t, θ) ≤ 0 (A.7d)
With initial conditions:
x(t0 ) = x0 (A.7e)
where φ is the objective function, x and y are time dependent differential and algebraic variables, u are the
control variables and θ are the time invariant parameters. Equation (A.7b) describes the DAE system
that the objective function, equation (A.7a), is subjected to with initial conditions of the differential
equations described in equation (A.7e). Equation (A.7c) and (A.7d) describes the equality and inequality
constraints of the system.
A.2.2 Dynamic Optimization and Discretization Approaches

There are two main classes of solution approaches for dynamic optimization problems: sequential methods
and simultaneous methods. Sequential methods use control vector parametrization to make the problem
tractable [59]. The (variable length) time horizon is divided into several control intervals and the control
variables are changed according to a predetermined function. Examples of control profiles can be seen in
Figure A.1 and include polynomials with continuity, piecewise linear profiles with and without continuity
constraints, and piecewise constant. Sequential methods are easy to set up and ensure accuracy due
to robust variable step integration algorithms. Nevertheless, a time consuming step for integrating
sensitivities of the objective function and constraints is required. In simultaneous optimization methods,
the state and control profiles are both discretized a priori (e.g., by collocation of finite elements). This
method leads to large non-linear programming problems and needs efficient algorithms to be solved [60].
51
Polynomial with
continuity
Piecewise
linear with
continuity
u(t)
Piecewise
linear
Piecewise
constant
t
Control Interval
Figure A.1: Examples of control profiles, adapted from [38]
A.2.3 Solution Algorithms

Many different solution algorithms have been implemented to solve nonlinear optimization problems. One
of the most frequently used is the sequential quadratic programming method. The objective function is
approximated with a quadratic equation, and the constraints are typically linearized. A search direction
is then determined which seeks to minimize the objective function. This method has proven effective
for both small and large nonlinear optimization problems. However, the method only guarantees locally
optimal solutions [61].
A.3 Reduced Order Modeling for Scheduling Applications

Plant schedules determine the order, rate, or quality of the products made in a plant. The schedule
is optimized to minimize off-spec product and inventory costs. For example, a chemical batch reactor
can produce several products depending on the operating conditions in the reactor or the mix of inputs
into the reactor. The optimal scheduling problem then seeks the manufacturing order which minimizes
production time and maximizes profit while satisfying demand for each product [2].
A.3.1 Scheduling
An optimal schedule can be determined using a static or dynamic plant model. Conventional approaches
use a static plant model. The static model incorporates the dynamics of the process by including a table
of transition times, τs , which are determined off-line. Allowable transitions between operating states
in the process are also included in the table. The table guarantees that no constraints are violated.
However, it is beneficial from a profitability perspective to account for the detailed process dynamics if
multiple products are produced in order to determine the minimal (feasible) transition time. This calls
for a close integration of scheduling with the control system to maximize profit and minimize the effect
of disturbances on the system [6].
Full dynamic scheduling uses a detailed dynamic plant model in the optimal scheduling problem.
A full dynamic scheduling approach can describe the dynamic behavior during transitions and likely
find the optimal production schedules, while ensuring that operational constraints are met. However,
the optimization problem can quickly become exceedingly large and difficult to find a solution in a
52
reasonable time frame. To overcome this problem, a reduced order model can be used to represent the
dynamic behavior of important process variables [2, 6].
A.3.2 Reduced Order Modeling

Park et al. [8] obtained promising results when comparing the optimal schedule determined using a full
dynamic model and a reduced order dynamic model. The reduced order model performed as well as the
full dynamic model when the plant model accurately matches the actual system, and outperforms the
full model when plant-model mismatch exists. In the work by Park et al. [8], the low order model was
established analytically based on the input-output behavior of the feedback control system of the plant.
If an input-output relationship cannot be established analytically, system identification can be used to
find a data-driven reduced order model [2, 6].
A.3.3 Model Identification

System identification consist of finding an input-output relationship from observed (historical) data
and/or physical knowledge about the system. These models can be “white box”, “black box”, or “gray
box”. In a “white box” model, the system is completely known and can be described by first principles
physics. A “black box” model, on the other hand, contains no physical insight of the system and is
solely based on observed data. Examples of linear “black box” models are auto regressive models and
transfer function models. Nonlinear “black box” models include artificial neural networks and nonlinear
Hammerstein-Wiener models. “Gray box” models contain either some physical knowledge about the
system combined with a “black box” model or physical insight into the system is used to choose a certain
kind of “black box” model [62].
Four basic steps needs to be done to find an accurate “black box” model: i) establish a set of historical
data, ii) find a set of candidate models, iii) have a criterion of fit, and iv) establish validation data and
validate the model [63]. A simple example of a linear “black box” model is the transfer function model.
A transfer function is the ratio between the Laplace transformed output and input of a system. The
Laplace transformation of x with respect to t is defined as [4, 64]:
Z ∞
X̄ = xe−st dt (A.8)
0
where s is a complex variable. The advantage of using a Laplace transformation is that a linear differential
equation problem is transformed into an algebraic problem. It can then be manipulated using algebraic
laws. A transfer function is described as:
Y (s)
G(s) = (A.9)
X(s)
where Y (s) and X(s) are usually described by polynomials in the Laplace domain. The roots of X(s) = 0
and Y (s) = 0 are called the poles and zeros of the transfer function, respectively. The number of
poles must be greater then the number of zeros for the transfer function to be physically realizable.
Furthermore, the number of poles describes the order of the system [4, 65].
A transfer function is used to describe the dynamic behavior of a system and relates one input, x, to
one output, y. The order of the transfer function corresponds to the order of the differential equation
which it represents. A series of transfer functions is commonly used in process control when relating the
input signal to the output through several “blocks” or sub-systems (e.g., relating the valve position to the
flow rate, and the flow rate to the outlet temperature), see Figure A.2. The resulting transfer function
is then the product of the two transfer functions [4]:
G(s) = G1 ∗ G2 (A.10a)
y = G(s)x1 (A.10b)
53
x1 G1 x2 G2 y
Figure A.2: Transfer functions in series
Transfer functions in parallel (e.g., Figure A.3) can be added together to obtain the overall transfer
function:
G(s) = G1 + G2 (A.11a)
y = G(s)x (A.11b)
G1
+
u y
+
G2
Figure A.3: Transfer function in parallel
An example of a non-linear “black box” model is the Hammerstein-Wiener model [66]. Such models
consists of a linear transfer function dynamic block with static nonlinear blocks which transform the
inputs and outputs, as seen in Figure A.4. The nonlinear blocks can have different structures, e.g.,
piecewise linear, networks, or polynomial. An example of a Hammerstein-Wiener model with a second
order polynomial as the input and output nonlinear blocks is presented below:
u0 = Au2 + Bu + C (A.12a)
0 0
y = Gu (A.12b)
y = ay 02 + by 0 + c (A.12c)
where u is the input to the model, y is the output, G is the transfer function model, and A, B, C, a, b,
c are parameters.
u u’ y’ y
Input Non- Output Non-
Linear Block
Linearity Linearity
Figure A.4: Hammerstein-Wiener model
54
Appendix B
Results
B.1 Steady state optimization

The optimal values of the decision variables at every production rate, see Table 4.1, are presented in
Figure B.1. Polynomial curves were fitted to this data for the purpose of calculating the optimal values
of the manipulated variables at continuously variable production rates. A simple linear model was fitted
to the feed flow data, F1,4 , see Figure B.1c. A second order polynomial was fitted to the fraction of
product refluxed to the condenser, rgasdraw , and the fraction of air liquefied in the PHX, kwithdrawn ,
with an R2 value of 0.97 and 0.98, respectively (Figure B.1b and B.1a).
0.5117 0.032
Data points Data points
Fitted polynomial 0.03
Fitted polynomial
0.5116
kwithdrawn [mol N2 reflux/mol N2 total]
0.028
rgasdraw [mol liquid air/mol air total]
0.5115
0.5114 0.026
0.5113 0.024
0.5112 0.022
0.5111 0.02
0.511 0.018
0.5109 0.016
0.5108 0.014
16 17 18 19 20 21 22 23 24 16 17 18 19 20 21 22 23 24
Production rate [mol/s] Production rate [mol/s]
(a) Fraction liquid air entering the column (b) Fraction of product entering the condenser
48
Data points
46
Fitted polynomial
44
42
F1,4 [mol/s]
40
38
36
34
32
30
16 17 18 19 20 21 22 23 24
Production rate [mol/s]
(c) Feed Flow
Figure B.1: Optimal control variables for different production rates with polynomial approximations
B.2 Outputs From the Scheduling Using the Reduced Order Model
In Figure B.2 the U A value, the reboiler liquid level and the production rate expected from the scheduling
using the reduced order model is given for 72 hours.
55
3 1.5
2 1
Reboiler Liquid Level

1 0.5
UA
0 0
−1 −0.5
−2 −1
−3 −1.5
0 20 40 60 80 0 20 40 60 80
(a) UA (b) Impurity

1.5
0.5
Production
−0.5
−1
−1.5
0 20 40 60 80
Time [hour]
(c) Production
Figure B.2: Expected variable and constraint outputs from the scheduling using the reduced order model for 72
hours
56
www.kth.se

Air Separation Unit

Uploaded by

Copyright:

Available Formats

Air Separation Unit

Uploaded by

Document Information

Copyright

Available Formats

Share this document

Share or Embed Document

Sharing Options

Did you find this document useful?

Is this content inappropriate?

Copyright:

Available Formats

Air Separation Unit

Uploaded by

Copyright:

Available Formats

DEGREE PROJECT, IN CHEMICAL ENGINEERING , SECOND LEVEL

STOCKHOLM, SWEDEN 2015

Integrated Scheduling and Control of

KTH ROYAL INSTITUTE OF TECHNOLOGY

SCHOOL OF CHEMICAL SCIENCE AND ENGINEERING

A Master’s Thesis at KTH Royal Institute of Technology

Supervisor: Dr. Michael Baldea*

July 13, 2015

*University of Texas at Austin, McKetta Department of Chemical Engineering

**KTH Royal Institute of Technology, Department of Chemical Engineering and Technology

3 Full Dynamic Model of an Air Separation Unit 7

4 Integrated Scheduling and Control Using Reduced Order Models 26

5 Results and Discussion 35

6 Conclusions and Future Work 44

1.1 The hierarchy of decision making time scales in a chemical plant . . . . . . . . . . . . . . 2

2.1 Flow sheet of the process . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6

3.1 Overview of the column . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8

4.1 Scheduling procedure . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26

A.1 Examples of control profiles . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52

3.1 Indices, variables and parameters used in the column model . . . . . . . . . . . . . . . . . 11

4.1 Production levels used in the steady state optimization . . . . . . . . . . . . . . . . . . . . 27

2.1 Air Separation Unit

2.1.1 Electricity Pricing Options for ASUs

2.1.3 Dynamic Operation and Scheduling of ASUs Subject to Time-Varying

2.1.4 Control Strategies of ASUs

2.1.5 Dynamic models of ASUs

Primary Heat Exchanger

Figure 2.1: Flow sheet of the process

Full Dynamic Model of an Air

3.1 Flow sheet overview

3.2 Limitations and Assumptions

Vi , yi,j Li-1 , xi-1,j

Vi+1 , yi+1, j Li , xi,j

V30 , y30,j L29 , x29,j

Figure 3.1: Overview of the column

Overall mass balance

Component mass balance

enthalpies were computed as:

dTi dxi,j ∂hL ∂hL

Now equation (3.11) is an algebraic expression containing Vi .

3.3.2 Integrated Reboiler and Condenser

Figure 3.2: Overview of the integrated reboiler/condenser

Qcondenser = U A(Tcondenser − Treboiler ) (3.19)

where pcondenser is the pressure of the condenser [bar] and psat

where sp denotes set point and Kc is the controller gain.

Component Mass Balance

dxreboiler,j L30 (x30,j − xreboiler,j ) − Vreboiler (yreboiler,j − xreboiler,j )

3.3.3 Primary Heat Exchanger

Segment 1 Segment n Segment 50 Segment 1

Air Feed Liquid air feed

Vapor air feed

Figure 3.3: Sideview of the PHX with the turbine

F1,1 , hV1,1,n F1,1 , hV1,1,n+1

F1,2 , hV1,2,n F1,2 , hV1,2,n+1

F1,3 , hV1,3,n F1,3 , hV1,3,n+1

F1,4 , hV1,4,n-1 F1,4 , hV1,4,n

Figure 3.4: One segment in zone 1

Metal wall temperature

Figure 3.5: Overview of the PHX

Vjlayer = (W − 2tedge )Nj Vjslot (3.30)