Chapter Overview: 2: Modeling Basics

CHAPTER OVERVIEW
2: MODELING BASICS
2.1: Verbal Modeling- process description, control specifications, and connections
2.2: Degrees of Freedom- importance, calculation procedure, and examples
2.3: Incidence Graphs -Interpretations, Consistency, and Inconsistency
2.4: Excel Modeling - logical models, optimization with solver for nonlinear regression, sampling random numbers
2.5: Noise Modeling - White, Pink, and Brown Noise, Pops and Crackles
2.6: Numerical ODE solving in Excel- Euler’s method, Runge Kutta, Dead time in ODE solving
2.7: Solving ODEs with Mathematica- How to find numerical and analytical solutions to ODEs with Mathematica
2.8: Fitting ODE parameters to data using Excel- Using regression to fit complex models in Excel
2.9: Helpful Mathematica Syntax- Hints on how to use Mathematica to model chemical processes
This page titled 2: Modeling Basics is shared under a CC BY 3.0 license and was authored, remixed, and/or curated by Peter Woolf et al. via source content
that was edited to the style and standards of the LibreTexts platform; a detailed edit history is available upon request.
1
2.1: VERBAL MODELING- PROCESS DESCRIPTION, CONTROL
SPECIFICATIONS, AND CONNECTIONS
Authors: (September 8, 2006) Brian McQuillan, Crystal Miranda, Brandon Quigley, and John Zhang
Stewards: (September 5, 2007) Kevin Luchi, Mike Nappo, Matt Neff, Lisa Schauman
2.1.1: INTRODUCTION
Every process requires a great deal of planning in order to successfully accomplish the goals laid out by its designers and operators. In order
to accomplish these goals, however, personnel that are not familiar with the design must fully understand the process and the functions of
the control systems. Control systems consist of equipment (measuring devices, valves, etc.) and human intervention (plant operators and
designers). Control systems are used to satisfy three basic needs of every process:
1. Reduce the influence of external disturbances
2. Promote the stability of the process
3. Enhance the performance of the process
Verbal modeling is used for creating and also understanding a process control system. Verbal modeling consists of first receiving and then
gathering information about the process. A step-by-step process is then used to describe the control systems used to satisfy constraints and
objectives that have been outlined. In the following sections you will read what requirements are generally outlined for process control and
the step-by-step method used to meet these requirements.
2.1.2: PREREQUISITE INFORMATION REGARDING A PROCESS

For the sake of this article it is assumed that a process has already been designed and that certain restraints and criteria are provided by
either a customer, management, or the government. The goal of this section is to classify the types of criteria that are usually given. These
criteria then become the conditions that the control systems employed must satisfy. In general, there will be five sets of criteria, often
coming from different people and institutions. By gathering all of these criteria you will be able to describe the control system. If you do not
have a complete list of these criteria you must research the process to determine these constraints before beginning the step-by-step process
below.
Safety
The safe operation of a process is the biggest concern of those working in the plant and those that live in the surrounding community. The
temperatures, pressures, and concentrations within the system should all fall within acceptable limits, and these limits can be dictated by
either government agencies or company policy.
Production Objectives
The production objectives usually include both the amount and purity of the desired product. This criterion is generally set by the company
or customer.
Environmental Regulations
These come in the form of restrictions on the temperature, concentration of chemicals, and flowrate of streams exiting a plant. State and
federal laws, for instance, may dictate the exit temperature of a cooling water stream into a lake in order to prevent harm to aquatic wildlife.
Operational Constraints
Equipment found in the plant may have their own unique limitations, such as temperature or pressure that require proper control and
monitoring. For instance, a thermocouple may be damaged at extremely high temperatures, thus the location of the thermocouple must be
accounted for.
Economics
In general, a company will operate so that its profits are maximized. The process conditions that maximize these profits are determined by
way of optimization. Many costs must be considered when optimizing process conditions. Some of these costs are fixed, or will not change
with process variables (i.e. equipment costs) and others are variable, or do depend on process variables (i.e. energy costs). The overall
process is usually limited by certain factors including availability of raw materials and market demand for the final product. Therefore, the
economics of a process must be well understood before process changes are enforced.
2.1.3: STEP-BY-STEP METHOD FOR DESCRIBING CONTROLS AND THEIR PURPOSE
2.1.1 https://eng.libretexts.org/@go/page/22542
1.) Describe the Process
A brief description of the general process is needed while not dwelling on the details and calculations involved. The major steps of the
process, as well as inputs and outputs of the process, should be stated. A simple diagram should be provided detailing the chemical process
to help visualize the process.
2.) Identify Process Objectives and Constraints
The objectives and constraints of the process must be identified before process control actions can be performed.
The process objectives include the type, quantity, and quality of the product that is to be produced from the process. The economic
objectives, such as the desired levels of raw material usage, costs of energy, costs of reactants, and price of products, should also be
identified.
The process constraints include three different categories: operational, safety, and environmental limitations. Operational constraints refer
to the limits of the equipment used in the process. For instance, a liquid storage tank can only hold a certain volume. Safety constraints
describe the limits when the people or the equipment may be in danger. An example would be a pressure limitation on a reactor, which if
exceeded, could result in an explosion. Environmental constraints limit how the process can affect the immediate surroundings. For
example the amount of harmful chemicals that can be released before damage is done to nearby water supplies. All of these constraints
should be mentioned to build a robust control system.
Careful reading of the information provided to you by the customer, management, and government is required in order to properly identify
each constraint and objective. Often times, the process objectives will be very clearly laid out by the needs of the customer or management.
Operational constraints, or the limitations of the equipment being used, must be researched for each piece of equipment used in the process.
Generally, by satisfying the operational constraints a good portion of safety constraints are satisfied as well, but additional safety constraints
may exist and must be investigated by researching company policy and governmental regulations. Environmental regulations also have to
be researched through resources such as the EPA and Clean Air Act. Satisfying the economic aspect is largely determined by manipulating
additional variables after all other constraints and objectives have been met.
3.) Identify Significant Disturbances
Disturbances, in the sense of process description, are defined as inputs or external conditions from the surrounding environment that have
certain properties that cannot be controlled by the plant personnel. Examples of disturbances include ambient air temperature, feed
temperature, feed flow rate, feed composition, steam pressure changes, and cooling water temperature changes. Disturbances can drastically
affect the operation of a unit. A control system should be able to effectively handle all process disturbances. As such, all possible
disturbances must be identified and these disturbances need to be accounted for by the development of contingency plans within the
process.
4.) Determine Type and Location of Sensors
A proper design must ensure that adequate measurements of the system are obtained to monitor the process. To meet this goal, sensors must
be chosen to accurately, reliably, and promptly measure system parameters. Such parameters include temperature, flow rate, composition,
and pressure. Placement of sensors is important both in the usefulness of measurements as well as the cost of the system. Sensors should be
placed such that the measured quantities are appropriate in addressing control objectives.
5.) Determine the Location of Control Valves
Valves must be placed in a location to control variables that impact the control objectives. For example, control of the temperature of a
reactor could be obtained by placing a valve on either the stream of heating / cooling fluids or by placing a valve on the feed stream to the
reactor. One must determine which streams should be manipulated to meet process objectives.
6.) Perform a Degree of Freedom Analysis
The degrees of freedom in a system are equal to the number of manipulated streams (determined in step 5) minus the number of control
objectives and control restraints (determined in step 2). A degree of freedom analysis is used to determine if a system is being under- or
over-specified by the process objectives. The degrees of freedom come from the number of knowns and unknowns that are specified within
the system. If there are extra degrees of freedom present in a system, unused manipulated variables can be used to optimize the process. If
there are negative degrees of freedom, a system is over-specified because more objectives and restraints exist than manipulated streams. In
this case, all objectives cannot necessarily be met simultaneously and the least important objectives must be neglected. A system with zero
degrees of freedom is fully specified. All objectives can be met, but there is no room for optimization.
7.) Energy Management
In any system with exothermic or endothermic reactions, distillation columns, or heat exchangers, energy management becomes a factor that
must be accounted for. Heat must be removed from exothermic reactions in order to prevent reactor runaway, and heat must be supplied to
endothermic reactions to ensure desired production rates. Strategies such as pre-heating feed streams with the excess heat from a product
stream are helpful in maintaining efficient usage of energy, however, they also result in more complex processes that may require more
intricate control systems.
8.) Control Process Production Rate and Other Operating Parameters
The production rate can be controlled by a variety of manipulated variables. One manipulated variable may be the feed rate. The plant feed
rate can be changed and each subsequent unit can use its controls to accommodate this change, ultimately resulting in a change in the final
production rate. Other manipulated variables may also include reactor conditions, such as temperature and pressure. Temperature and
pressure affect reaction rates and can be used to alter the final production rate. It is important to choose the most suitable manipulated
variable to control production rate.
In addition to the production rate, other control objectives must be effectively managed by manipulated variables. For example, temperature
of an exothermic reactor may be controlled by the flow of a coolant stream passing over it in order to avoid dangerous high temperatures.
The pressure of a reactor may be controlled by the flow of feed gas in order to comply with the pressure limitations of the vessel.
9.) Handle Disturbances and Process Constraints
The effects of disturbances should be minimized as much as possible, in order to maintain the system at desired conditions and meet all
process objectives and constraints. Feedback or feedforward are specific control techniques and are common ways to overcome
disturbances. A feedback control works by studying the downstream data and then altering the upstream process. The actions executed are
reactive. Feedback can be viewed as an if-then statement: if a feed's temperature is detected to be lower than desired, then steam can be used
to preheat the feed. Feedforward is a more proactive approach in that it adjusts a manipulated variable before the disturbance is felt in the
process. Hence, if a sensor indicates low temperatures upstream of the feed, the feedforward control will counteract the effect of the cooler
upstream temperatures by preheating the feed before the feed temperature is effected. Note that a disturbance must be detectable and
measurable in order for the feedforward control to fix the anticipated disturbance before the system is effected.
Additionally, if constraints are reached during the process, controls should be implemented to avoid safety, operational, or environmental
hazards. This can also be done with feedback and feedforward controls on manipulated variables.
10.) Monitor Component Balances
Every component within a process, whether it is inert or not, should be accounted for at every step of the system in order to prevent
accumulation. This step is more crucial in processes that involve recycle streams. If such a stream is present, a purge stream is often
necessary to remove unwanted components. In addition, component balances are used to monitor yield and conversion or reveal locations in
the process where loss may be occurring. In order to monitor component balances, composition sensors are used.
11.) Control Individual Unit Operations
Most systems used today in industry employ the use of multiple unit operations. Each of these unit operations, however, needs to be fully
controllable in the sense that it has a control system that can adjust manipulated variables in order to maintain other parameters. For
instance, if an absorber is present, the system must be able to control the liquid solvent feed as some ratio to the gas feed. Another example
is a crystallizer. The refrigeration load of the crystallizer must be controllable in order to control the temperature.
12.) Optimize the Process
In most cases, there will be certain aspects of a process that will not be dictated to a designer and can be changed to make the overall
process more economical for the company. These are referred to as "unaccounted for" degrees of freedom and can be implemented as new
control valves or adjustable controller setpoints.
2.1.4: ALTERNATIVE METHOD OF VERBAL MODELING

This alternative method is also described in Peter Woolf's Recorded_Lectures Lecture 3. For purposes of this description, it is scaled down
to a single unit process. However, this method can easily be applied to describe an entire system of unit processes.
1. Describe the process in words
Some of the important questions to answer before delving deeper into a model are:
What are the components entering the system?
How do they enter? Separately? Combined stream? What physical states are they in?
What happens inside the unit process and what comes out at each exit point?
Remember to keep this part simple. There is no need to include chemical formulations or equations of any sort. Just lay out the basic flow of
material.
2. Define the primary goal of the process
The primary goal should be simple. Often, it is to maintain a specific measured variable above a minimum or below a maximum. In this
step, the only thing that needs to be determined is what the main goal is, and a few supporting details about why this is an important goal to
achieve.
For example, a primary goal could be to minimize the concentration of Compound Y in orange juice because studies show Compound Y
gives the juice a bad aftertaste.
3. Identify secondary processes that influence the primary goal
In a typical unit process, the primary goal will be directly influenced by one or two other aspects of the system. These can include
temperature, pressure, inlet conditions, and more and can occur at various points in the process.
The goal of this step is to determine which of these other process variables will be most likely to influence the primary goal and to step
down from there.
For example, the temperature of the orange juice mixer could have the greatest influence on production of Compound Y.
4. Identify safety and environmental risks
Next, you need to identify all of the points in the process that represent any type of risk. This will be important later in determining which
system variables need to be monitored.
Step through your process and identify any points that pose a significant risk of the hazards shown in the following figure.
Examples include: Boilers represent fire and explosion risks. Any stream with a dangerous chemical can represent corrosive, poison,
environmental, or all three risks.
5. Identify major costs associated with the process
How much something costs to produce is obviously a big deal in manufacturing. Identifying the largest sources of cost is critical in finding
ways to reduce cost overall. Typical places to start identifying costs are at inlet streams (what is the cost of raw materials) and at any portion
of the process where heat is added or removed.
It is important to include the high costs that can be associated with the risks identified in Step 4. Often the high cost of failure and risk
exposure will determine what other seemingly costly steps must be taken to ensure the safety of the process.
6. Identify variables you can directly manipulate
The basics of the process have been laid out, and now it's important to determine what variables you can actually control. Typically, you
only have direct control over the simplest of variables: switches and valves. Essentially, this just means that you cannot, in fact, choose a
temperature for your system and implement it. What you can do, is control a valve or switch that activates heating or cooling to control the
temperature.
During this step, you must decide where it is important to place these valves and/or switches. Use the information acquired previously about
the primary goal and secondary effects to determine what variables are worth controlling. Remember that you don't need to put valves
EVERYWHERE! Control valves are not costless and can also add unwanted complexity to your system. If you need to isolate equipment,
you can install manual valves. Keep the control valves to the needed level.
7. Identify sources of variation
In order to write a control scheme, you need to know what values in your system will change and why. Some common causes of variation
include:
Environment: ambient temperature
Other processes upstream or downstream: variable inlet conditions or outlet demand
Economic forces: product worth, material costs
Operators
Identifying what aspects of your process can be affected by these forces will allow you to assemble a more complete control scheme.
8. Describe your control system in words
Before you start trying to write everything out in computer code and mathematical equations, take the time to lay out your controls in
words. This is much like preparing an outline before writing a paper. It can save you from many headaches later on.
One example of generic, simple syntax for verbal modeling is: Maintain [system variable] at specified level by adjusting [variable I can
control].
2.1.5: THE BARKEL METHOD OF VERBAL MODELING

This method is an elaboration upon the steps outlined in Mr. Barry Barkel's lecture on September 29, 2009.
"If you design a system, you have the ethical responsibility to control it."
1) Understand the Process
Before you can control anything, you have to understand the process and how different parts of the system interact with each other. Make
sure that the overall process is understood. This includes inputs and outputs as well as major steps, however specifics are not necessary. You
should also be able to construct a diagram to help explain the process.
2) Identify Operating Parameters
Operating parameters can include temperature, pressure, flow, level, etc. Choose the parameters to manipulate in your system that will
safely result in the desired output.
3) Identify hazardous conditions
Consider all possible dangerous aspects of your process when designing your system. This could include a chemical overheating or a vessel
overflowing. It is imperative that you ensure the safety of your operators and being aware of all hazardous conditions can aid in this.
4) Identify measurables
The main three measureables addressed in this class are: temperature, pressure, and flow. However, there are many more measureables,
some more common than others. Here are some more: pH, humidity, level, concentration, viscosity, conductivity, turbidity, redox/potential,
electrical behavior, and flammability.
5) Identify points of measurement
It is important to place sensors in locations so that efficient and reliable measurements are taken of the system to monitor the process. For
example, in a distillation column, temperature sensors will display different temperatures at different locations down the tower. The sensors
must be positioned so that accurate readings of the system are given. Also, it is necessary to place sensors in an area of constant phase.
6) Measurement methods (thermo couple? choose for range)
After identifying what is to be measured and where it will be measured, you have to decide how it will be measured. For example, a
thermocouple can be used to measure the temperature. When choosing the equipment, be sure to check that the conditions of use fall within
the recommended range of operation for the equipment.
7) Select control methods
Decide whether to use feedback, feed-forward, cascade, or other types of control methods.
8) Select control system
A control system is a set of devices that will manipulate the actions of other devices in the system. A control system ranges from having an
operator manually open and close a valve to running a system with feedback such as with PID controllers. Control systems can vary from
relatively cheap to expensive. When picking a system, it would be most economic to choose the cheapest one that gets the job done. An
example of a control system is a PIC, or programmable interface controller.
9) Select control limits
When choosing setpoints for controllers, you will also have to decide on a range that the values are allowed to fluctuate between before a
corrective change is made. When selecting these limits, keep in mind that "equal" and "zero" do not exist due to the infinite number of
decimals that electronics are now able to handle. So you must further define what "equal" means, when is it "close enough" to count as
"equal" And when is the number small enough to count as "zero", is 0.1 or 0.01 or 0.0000000001 count as "zero"?
10) Define control logic
As every process is different, a customized code for each process must be written to tell the system what to do. For example, when a level
control is a tank has reached a critically high point, the logic should spell out the necessary changes needed to bring the tank level back
down. For example, this could be partially closing a valve upstream of the tank or partially opening a valve downstream of the tank.
11) Create redundancy system
In the real world, you must balance cost and efficiency/safety. On one hand, you don't want an out-of-control system if one control fails. But
on the other hand, you can't afford to order two of everything. The critical point to keep in mind is to optimize the safety while minimizing
the cost.
12) Define "fail-safe"
A fail safe is a set up in the control logic to ensure that in the event of a failure of a control method, the system will automatically reach a
safe condition so that there is little to no harm done to other equipment or personnel.
13) Set lead/lag criteria
Valves and other equipment do not necessarily open/close or turn on/off at the exact instant a button is pressed or the control logic kicks in.
There is often lag time associated with each controller. You must determine how long this lag time is so that you can take it into account.
14) Investigate effects of change before/after
Be sure to investigate effects of changing each controller. For example, what are the effects of closing/opening this valve?
15) Integrate all systems
Ensure that all systems are working together and that there are no holes in the system. Make sure that information does not fall through any
cracks in the system.
[Note - we can use an example of the Barkel method - RZ]
2.1.6: COMMON ERRORS

1. Impossible direct manipulations e.g., Change the concentration of salt in a tank
2. Missing the forest for the trees e.g., Sacrificing product quality for tight level control on a tank
3. Excessive or insufficient control e.g., Control every variable because you can or ignore the possibility of significant disturbances
2.1.7: WORKED OUT EXAMPLE 1

PLEASE NOTE: ALL VALUES IN THE FOLLOWING PROBLEMS ARE FICTICIOUS BUT MEANT TO BE LOGICAL
A heat exchanger uses steam to heat a stream of water from 50ºF to 80ºF. The water enters at a flow rate of 20 gallons per minute
from a nearby lake. The process costs $65 per hour and yields a profit of $2 per gallon of product. The steam is provided by the
plant and is 1000ºF, which is a temperature that the pipes can sustain. For safety reasons, the exchanger may only run for 12
consecutive hours and requires 4 hours to cool down. Using more than 10,000 gallons of water per hour would cause an
environmental disturbance to the water source. The diagram is shown below. Verbally model this system.
== Solution ==
The purpose of the process is to heat an incoming stream of water from a temperature of 50ºF to a temperature of 80ºF. The main equipment
involved is a shell-and-tube heat exchanger.
The product specification of the process is water at a flow of 20 gallons per minute and a temperature of 80ºF.
Economically, the process costs $65 per hour to operate. There are no costs for the raw materials, as the only inputs to the system are water
and steam. The finished product produces a profit of $2 per gallon. The economic objective is to reduce process costs while producing
sufficient product.
The operational constraints and safety concerns are due to the pipes. The pipes can only sustain a temperature of 1000ºF. Safety is a concern
because attempting to heat the incoming water to a certain temperature may cause the heat exchanger to malfunction, leading to equipment
damage and possible burn injuries to nearby personnel. The system may only operate for 12 consecutive hours, after which the system will
need to be cooled down for 4 hours to avoid the aforementioned hazards. A simplified assumption is that there are no constraints on steam
because it is provided by the plant and causes no safety issues. The only environmental constraints involve the incoming water stream. The
incoming water is gathered from the nearby lake, and a stream of greater than 10000 gallons per hour would cause a disturbance in the
equilibrium of the lake.
Significant disturbances can be found in the ambient air temperature, variable flow rates of the feed, and the temperature of the steam.
4.) Determine the Type and Location of Sensors
A flow sensor (FM) is placed at the incoming water stream. A temperature sensor (TS) is located on the product water stream. A flow
sensor is not needed for the steam stream for this problem because this value is not needed for control. A sensor could be placed here but the
information is not needed for this problem.
A flow valve is placed at the entrance of the incoming water stream. A flow valve is placed at the entrance of the steam.
6.) Perform a Degree-of-Freedom Analysis
There are two manipulated variables: the flow of the water feed stream and the flow of the incoming steam. There are two control
objectives: the flow of the feed stream, monitored by the flow sensor, and the temperature of the product, monitored by the temperature
sensor. Therefore the system has zero degrees of freedom.
The incoming steam is used to transfer heat to the cool water feed. The temperature sensor on the product stream determines the applicable
setting on the steam flow valve.
The process production rate is controlled by the flow valve on the entering water stream. The water temperature is controlled by the flow
valve on the incoming steam.
Changes in the ambient air temperature can be detected by the temperature sensor, and can be corrected by the flow valve on the incoming
steam stream. Variable flow rates of the water feed stream can be detected by the flow sensor and compensated by adjustments on the flow
valve on the water feed stream. Changes in the temperature of the steam can be detected by the temperature sensor. The flow valve on the
steam stream can be adjusted to increase or decrease the flow of steam and subsequently the amount of heat exchanged.
A vent is located on the heat exchanger to release excess steam from the system. Aside from that, any accumulation is unlikely and can be
neglected.
The outlet temperature of the product stream is controlled by the flow valve on the steam feed stream. The flow of the incoming water
stream is controlled by the flow valve on the water feed stream.
One might notice that the process is only using 1,200 gal/hr of water, well below the 10,000 gal/hr environmental constraint. If the profit of
the process is linear with the flow-rate of water, then increasing the flow-rate of water will increase the profits for the company. (With the
constraints specified, this is a Linear Programming optimization problem. The optimal setpoint falls on a boundary condition.) However, the
flow-rate of water entering the system is already specified, which results in zero degrees of freedom. (Zero degrees of freedom implies there
are no further control valves or setpoints.) Further investigation should be conducted to determine the reason for the flow-rate specification.
When considering increasing the flow-rate of water into the system, one should also check that the other constraints are not violated.
A process converts phenol into salicylic acid through a series of two reactors. Phenol and NaOH are fed in the liquid phase into the
first reactor where it reacts with gaseous carbon dioxide that is pumped in. Assume constant fresh feed temperature and that the
feed flow rate is within operational constraints. Management has dictated that salicylic acid production must be 200 moles per hour.
Also, management would like the product stream to have a molar composition of 80% salicylic acid. Due to environmental
concerns, a maximum flow rate of 10000 gallons per hour of cold water can be used to cool the first reaction chamber. The valve
controlling the flow of cold water does not allow a flow rate in excess of 7500 gallons of water per hour. The salicylic acid product is
used in another process to produce aspirin, which has a market value of $10 per mole. The first reactor can be operated at pressures
up to 200 atm, while the second can be operated at pressures up to 10 atm. The first reaction is exothermic, while the second
reaction is assumed to generate negligible heat. A diagram of this process is shown below, as well as the reaction scheme. Verbally
model this system.
== Solution ==
The purpose of the process is to convert phenol into salicylic acid, the precursor for aspirin. First, the phenol is reacted with gaseous carbon
dioxide (CO2) and sodium hydroxide (NaOH) under high pressure and temperature in the first reactor. The product of this reaction is then
combined with sulfuric acid to form salicylic acid in the second reactor. The reaction scheme is shown above.
The process is expected to produce 200 moles per hour of salicylic acid. The product stream must contain at least 80% by moles of salicylic
acid. The equipment used in the process dictates the operational limitations. The first reactor vessel can be operated up to a pressure of 200
atm, while the second reactor vessel has a 10 atm upward pressure limit. As such, pressures in excess of these limits must be avoided. Since
the first reactor will generate a significant amount of heat, the heat must be removed to avoid damage to equipment and possible runaway
reactions. Therefore, a heat exchanger (in the form of a reactor jacket in this case) with cool water should be included to decrease the
temperature of the reactor. Economic concerns demand that phenol, sodium hydroxide, and sulfuric acid should not be used in extreme
excess. The costs of these materials and the energy costs required to process them affect the overall profitability, so these compounds should
not be wasted. Environmental regulations limit the use of water to cool the reactor at 10000 gallons per hour, however the valve constraints
limits the amount of water to only 7500 gallons per hour.
The amount of cold water available to cool the reactor can be considered a disturbance because it comes from a reservoir outside of our
control. The ambient temperature is also a disturbance. If it drastically increases, the amount of cold water needed to cool the reactor would
need to increase as well. Composition of the feed streams will be assumed to be constant in this example. Therefore, they are not considered
disturbances.
4.) Determine the Type and Location of Sensors
A temperature sensor (TS) and pressure sensor (P) are located on the stream exiting the first reactor vessel. A flow meter (FM) is located on
the product stream leaving the second reactor. A composition sensor (CS) will also be located on the product stream leaving the second
reactor. The pressure drop can be controlled through the decompressor and thus is a control.
Control valves are located on the feed stream containing the phenol and sodium hydroxide, the incoming cold water to the first heat
exchanger, and the sulfuric acid feed stream. There is also a pump located on the carbon dioxide stream that enters the reactor.
6.) Perform a Degree of Freedom Analysis
There are 3 valves, 1 pump, and 1 decompressor but 5 objectives. This results in zero degrees of freedom. The valve located on the sulfuric
acid feed stream is meant to meet the composition constraint placed on the product stream leaving the second reactor. The valve located on
the feed stream carrying the reactants is set to satisfy production requirements. The valve on the cold water stream is used to maintain
reactor temperature, which satisfies an operational constraint. The pump is to ensure the correct pressure is achieved in the reactor, also
satisfying an operational constraint. The decompresser is to maintain a pressure of less than 10 atm in the second reactor, thus satisfying
another operational constraint.
The heat from the exothermic reaction in the first reactor is transferred to the cold water stream. The hot water stream exiting the reactor
vessel jacket could be used to heat streams on other processes. The second reactor is assumed to generate negligible heat during the
reaction, thus any release of heat from the reactor will be considered safe to release into the environment surrounding the process.
The production rate is measured by the flow sensor on the product stream and this signals the control valve on the feed stream through a
feedback mechanism to change the production rate as necessary.
If the temperature sensor on the reactor exit stream exceeds a certain level due to a diminished cold water supply, the feed stream valve
would decrease the amount of reactants entering the reactor. The amount of feed would also be decreased if more than 7500 gallons per hour
of cooling water were needed, as this is an operational constraint. If the pressure gauge controlling the pump begins to read higher than
allowed pressures, the pump would decrease the flow of the carbon dioxide entering the reactor. Also, if the pressure gauge reads out a
pressure that will be too high for the second reactor, the decompresser will be allowed to disperse more pressure. If ambient air temperature
drastically increases, the temperature sensor would open the cold water valve allowing more cooling water to enter the reactor vessel jacket.
If the composition of the product stream falls below 80 mole percent of salicylic acid, then the valve controlling the sulfuric acid feed would
allow more sulfuric acid into the second reactor to increase the conversion of reactants.
The composition sensor and flow meter on the product stream leaving the second reactor will account for every species to ensure that there
is no accumulation or loss within the system.
The first reactor vessel's pressure is fully controlled by the pressure gauge and pump system and its temperature is fully controlled by the
temperature sensor which controls the reactant feed valve and the cool water valve. The second reactor's pressure is fully controlled by the
same pressure gauge and the decompresser system, and its temperature will be highly dependent on the amount of cooling water used to
cool the product exiting the first reactor.
Since there are no unaccounted degrees of freedom, there are no valves to adjust in order to optimize the process. It should be noted,
however, that if there was no constraint on the composition of the product stream, the sulfuric acid feed valve would have become an
unaccounted for degree of freedom. If this had been the case, the valve could be adjusted to maximize the profit of the process. In order to
maximize the profit, the benefits of having higher conversion and more product would have to be weighed against the increase costs of
using more sulfuric acid feed.
2.1.9: SAGE'S CORNER

An example of verbal modeling
http://www.youtube.com/watch?v=YTs_IWccnvw&feature=player_embedded
Powerpoint Slides from the example
2.1.10: REFERENCES
Luyben, William L., Tyreus, Bjorn D., and Luyben, Michael L. “Chapter 8: Eastman Process” in Plantwide Process Control, McGraw-Hill,
New York, pp. 251-272.
Luyben, Michael L., Tyreus, Bjorn D., and Luyben, William L., "Plantwide Control Design Procedure" in AIChE Journal Dec. 1997, Vol.
43, No. 12 pp. 3161-3174.
Riggs, James B. and Karim, M. Nazmul. “Chapter 17: Multiunit Controller Design” in Chemical and Bio-Process Control, Ferret
Publishing, pp. 491-504.
Stephanopoulos, George. Chemical Process Control: An Introduction to Theory and Practice, Prentice Hall, New Jersey, pp. 1-41.
This page titled 2.1: Verbal Modeling- process description, control specifications, and connections is shared under a CC BY 3.0 license and was authored,
remixed, and/or curated by Peter Woolf et al. via source content that was edited to the style and standards of the LibreTexts platform; a detailed edit history
is available upon request.
2.2: DEGREES OF FREEDOM- IMPORTANCE, CALCULATION PROCEDURE,
AND EXAMPLES
Authors: (13 December 2009) Jon Butler, Matthew J LaVelle
2.2.1: INTRODUCTION
In control engineering, a degree of freedom analysis is necessary to determine the regulatable variables within the chemical process. These
variables include descriptions of state such as pressure or temperature as well as compositions and flow rates of streams.
2.2.2: DEFINITION
The number of process variables over which the operator or designer may exert control. Specifically, control degrees of freedom include:
1. The number of process variables that may be manipulated once design specifications are set
2. The number of said manipulated variables used in control loops
3. The number of single-input, single-output control loops
4. The number of regulated variables contained in control loops
The following procedure identifies potential variables for manipulation.
2.2.3: THE PROCESS

The method we will discuss is the Kwauk method, developed by Kwauk and refined by Smith. The general equation follows:
Degrees of freedom = unknowns - equations
Unknowns are associated with mass or energy streams and include pressure, temperature, or composition. If a unit had Ni inlet streams, No
outlets, and C components, then for design degrees of freedom, C+2 unknowns can be associated with each stream. This means that the
designer would be manipulating the temperature, pressure, and stream composition.
This sums to an equation of
Total Unknowns = Ni*(C+2) + No*(C+2)
If the process involves an energy stream there is one unknown associated with it, which is added to this value.
Equations may be of several different types, including mass or energy balances and equations of state such as the Ideal Gas Law.
After Degrees of Freedom are determined, the operator assigns controls. Carrying out a DOF analysis allows planning and understanding
of the chemical process and is useful in systems design.
2.2.4: APPLICATIONS
Single phase systems
All outlet streams have the same composition, and can be assumed to have the same temperature and pressure
Multiple phase systems
An additional (C-1) composition variable exists for each phase
Complete Process
When connecting units which share streams, one degree of freedom is lost from the total of the individual units
2.2.5: EXAMPLE 1: BLENDER

This example (from the ECOSSE Control Course web page listed below) investigates degrees of freedom in a simple vapor mixing unit.
Two gaseous streams enter a vessel and exit as a single well-mixed stream (Figure 1). We will apply the above equation to determine
degrees of freedom.
Figure 1: Blender Schematic
Here, there are 3 streams, each with C+2 unknowns for a total of 3C+6 Unknowns.
We have C mass balances and 1 energy balance for a total of 2C+1 equations. We also know composition, pressure, and temperature for the
incoming streams. This adds 2C+2 to the equation. Putting everything together gives:
Degrees of freedom = 3C+6 - (2C+1 + 2C+2)
Hence, the system has 3 degrees of freedom. Therefore, we can fix outlet composition, pressure, and flow rate. Figure 2 shows an example
control scheme:
Figure 2: Blender with controls
2.2.6: REFERENCES
Ponton JW, 1994, Degrees of Freedom Analysis in Process Control, Chemical Engineering Science, Vol. 49, No. 13, pp 1089 - 1095.
eweb.chemeng.ed.ac.uk/courses/control/restricted/course/third/course/formal.html
This page titled 2.2: Degrees of Freedom- importance, calculation procedure, and examples is shared under a CC BY 3.0 license and was authored,
remixed, and/or curated by Peter Woolf et al. via source content that was edited to the style and standards of the LibreTexts platform; a detailed edit history
is available upon request.
2.3: INCIDENCE GRAPHS -INTERPRETATIONS, CONSISTENCY, AND
INCONSISTENCY
3.1 Introduction
Incidence graphs are useful tools for chemical engineers to model various relationships in a process. They are used to systematically map an
entire chain of processes and controllers to describe the effect, which each element has on the others. This helps to visualize the possible
process pathways or a chain of effects.
Below is an example of an incidence graph. The circles are nodes which represent a particular device or controller, e.g. a temperature
sensor. The arrows indicate the directional pathway of influence. For instance, a temperature sensor node would be connected to an adjacent
node (such as a heat exchanger). This heat exchanger node would in turn be connected to other nodes representing devices or controllers.
This type of diagram could be extremely useful in identifying the redundancies within a control system. A more elaborate application of
incidence graphs will be discussed in the Worked Out Examples section.
2.3.1: MONOTONICITY
Before proceeding to the next few sections, it is imperative to understand the basics of monotone systems. A monotonic function defines a
function that preserves a given order. Monotonic functions are classified as monotonic increasing or monotonic decreasing. Taking the
example of a function f(x), for a monotonic increasing function an increase in variable ‘x’ would result in an increase in f(x). On the other
hand, for a monotonic decreasing function an increase in variable ‘x’ would result in a decrease in f(x).
Monotonic Increase : If x≥y, then f(x)≥f(y)
Monotonic Decrease : If x≥y, then f(x)≤f(y)

As long as a function is either always monotone increasing or monotone decreasing, it is referred to as a monotone function.
Monotonic functions are ideal systems to control because changes in one variable/device lead to known changes in other variables/devices.
A non-monotonic function on the other hand, does not yield a straightforward output for a particular input. In these cases IF-THEN-ELSE
statements are used to illustrate the controls for the nodes. For instances, IF node 1 is activated, THEN node 2 is repressed and so on.
An example of a possible non-monotonic graph
2.3.2: INCIDENCE GRAPHS

Incidence graphs (also “Levi” or “Causation” graphs) consist of monotonic connections between nodes, whereby each connection indicates
a particular effect of one factor on another. The term “monotonic” means that a change in the initiating device influences only one aspect of
the target device. This simplification facilitates describing how control systems work through the use of incidence diagrams.
In reality, not all interactions are monotonic. Some control devices can simultaneously affect several aspects of other devices or variables.
This may result in a non-monotonic relationship, in which one factor increases another over one region, but diminishes it over another
(imagine a parabola, which is a non-monotonic curve). Temperature often has such a relationship with the rate of a reaction: for an
exothermic reversible reaction, increasing the temperature will increase yield up to a certain point; however, further increase will diminish
the yield.
Incidence graphs can only be used to represent monotonic relationships. Therefore, over simplification of a system can cause problems in
constructing an incidence graph. Similarly a non-monotonic relationship can be resolved by breaking down a factor into separate
components, each of which has only monotonic relationships with other factors.
When two factors are connected by multiple signal paths, incoherence between the results of different paths leads to inconsistency. In other
words, the outcome cannot be predicted because the incidence graph does not communicate quantitative information. For instance, in the
reaction rate example, reaction rate could be separated into the variables which affect it: concentrations, rate constant, and equilibrium
constant. An incidence graph can then be constructed to demonstrate that temperature and reaction rate have an inconsistent relationship.
A diagram is consistent (to be discussed later in detail) if all signal paths between two nodes, for every pair of nodes, are in agreement with
one another (i.e. they have the same result).
2.3.2.1: TRAVERSING
Every incidence graph establishes a particular way to depict the direct effects of one device on its neighbors. Typically, a line will connect a
node to its neighbors--these lines are terminated by either a bar or an arrowhead.
An arrowhead at the end indicates that the initiator (from which the line is coming) is either increasing the value of the target (where the
arrow is pointing) or simply activating it (turning it on). The context of the connection will determine whether something is being increased
or turned on. For example, it would be illogical to assume that a drop in temperature would cause a temperature controller to deactivate a
pressure controller. Instead, it would be more prudent to assume that the temperature controller would adjust the pressure sensor.
arrowhead bar
A perpendicular bar indicates that the initiating device is either decreasing the value of the target device or inhibiting it (turn off). Note that
these indications are not necessarily strictly followed. Once again, context of the situation determines which interpretation should be used.
perpendicular bar
Graphs consisted of the arrowhead bar and the perpendicular bar allow you to see the result of an increase/decrease in one aspect of the
initial device on a designated device in the system. Everything done after the initial device comes as a result of changing the initial device.
In addition, these graphs provide the possibility of simplifying complicated relationships between devices into one conceptual unit. The
relationship between the initiating device and target nodes can be described by the following set --> {initiator,target1,target2,etc}. It is read
as "the increase or decrease in the initiating device affects the target nodes in the following manner... (+ for increase/activate and - for
decrease/inhibit). The initiator can be any device or variable in the system as long as the arrows support the designation. In other words, you
can show the effects of any device on any other device, as long as the pathway is supported by the diagram. There are often multiple routes
that a signal can travel.
2 helpful tips to keep in mind as you interpret an incidence graph:
1. Always travel from tail to head (or source to target) when traversing the graph from the initiating device to the designated device. Keep
in mind that it is possible to have arrowhead and/or perpendicular bars pointing in both directions in between two devices.
2. Arrowhead bars don’t always mean increase as perpendicular bars don’t always mean decrease. Arrowhead bars indicate that the effect of
performing an action to the device at the tail of the bar would yield the same effect to the device at the head of the bar. Penpendicular bars
indicate that the effect of performing an action to the device at the tail of the bar would yield the reverse effect to the device at the head of
the bar.
The following is a simple example of an incidence graph. We will refer to the circles below as “nodes” from this point onward. Nodes are
the selected devices or variables in a control system.
This incidence graph indicates that:

a. A change in an aspect of 1 will change some aspect of device 2 in the same direction -- otherwise written as {+,+}
The ends indicate the direction of influence as well as the type of influence.

a. increase in an aspect of 2 (initiator) decreases some aspect of 1 (target)
b. indicated by {+,-}
This next example builds on the previous two...

Given an arbitrary 3,1,2 order, an increase in 3 (initiator) will make 1 (target) go up. 1’s increase will then make 2 decrease. The symbolic
representation for this would be {+,+,-}. One can see that this diagram shows how node 3 influences everything directly and indirectly
connected to it.
Now, given a 1,2,3 order, node one becomes the initiator. The symbolic representation then changes to {+,-,-}. NOTE: the last node is
negatively affected because 2’s impact on 3 is shown to increase when 2 increases. But since 1 decreases 2, then the opposite is true!
The main objective of these diagrams is to determine the effect that a change in the starting node has on another node in the system.
However, in some cases, you can follow a pathway that ends at the starting node in order to check for ambiguity. For example, the 1,2,3,1
order yields a result of {+,-,-,-}. Since increasing node 1 cannot cause a decrease in itself, it is ambiguous.
Incidence graphs for complex controls systems will usually indicate that there are several ways for a device to influence another device in
the system.
Multiple paths
Remember these are done with respect to increases in a node’s value. So, if a path leads you to a point where you decrease a node, then you
have to do the opposite of what the diagram is indicating for all steps after that (for that route).
The incidence graph indicates:
1. For an order {1,2,3}: you observe {+,-,-}... remember not to count beginning node twice…
2. For an order {1,3}: you observe {+,+}
2.3.3: CONSISTENT, INCONSISTENT, AND PARTIALLY CONSISTENT

When designing controls for systems it is always important to keep in mind the need for redundacy and feedback control. Redundancy in a
system will ensure that a system is monitored by multiple controllers. So, if one controller was to fail there would be another to ensure the
same desired outcome. Consistent graphs can be used to describe redundancy. Feedback control is also very important and it is not good to
have the system overreact to external stimuli. A controlled response is desired. Inconsistent graphs can be used to describe feedback control.
Additionally, systems exist that both need to be redundant and feedback controlled and these can be described with partially consistent
graphs.
2.3.3.1: CONSISTENT GRAPHS

Consistent pathways happen when all pathways from the same starting node, that lead to the same target node, give the same output.
For example,
All possible pathways from node 1 to target node 4 yield the same output. These are called consistent pathways.
Paths(s) Sign series

1,4 (+,-)
1,2,4 (+,+,-)
1,3,4 (+,+,-)
From the table and graph above, all pathways from node 1 to node 4 yield the same outcome. Node 4 is being decreased.
A consistent graph is made up of only consistent pathways for all possible pathways between any combination of two nodes. The graph
below is an example of a consistent graph.
1,2,3 (+,+,-)
1,3 (+,-)
3,1,2 (+,-,-)
3,2 (+,-)
1,3,1 (+,-,+)
1,2,3,1 (+,+,-,+)
3,2,3 (+,-,+)
3,1,3 (+,-,+)
3,1,2,3 (+,-,-,+)
All pathways leading from node 1 to node 3 -

1,2,3 : (+,+,-)
1,3 : (+,-)
Since the two pathways cause the same change to node 3, this sub-pathway is consistent.
All pathways leading from node 3 and to node 2 -
3,1,2 : (+,-,-)
3,2 : (+,-)
Since the two pathways cause same change to node 2, this sub-pathway is consistent.
All pathways leading from node 1 and back to node 1 -
1,3,1 : (+,-,+)
1,2,3,1: (+,+,-,+)
Since the two pathways result in the same change to node 1, this sub-pathway is consistent.
3,2,3 : (+,-,+)
3,1,3 : (+,-,+)
3,1,2,3: (+,-,-,+)
Since the three pathways give the same effects to node 3, this sub-pathway is consistent.
Because all sub-pathways of this incidence graph are consistent, the entire process is consistent.
The table shown below gives examples of many different incidence diagrams and shows whether or not they are consistent.
Another type of consistent graph is one in which the valves of the system are not dependent on one another for feedback. In other words,
there is no proportion of the output signal from a particular valve of a system is passed to the input of another valve. An example of this
would be to consider the following process diagram:
This process is of a distillation column in which the bottom columns collect high intensity liquids and the liquid is removed from a series of
evaporators. The evaporator is a single-effect evaporator equipped with a simple heat exchanger, a separation vessel and a steam jet ejector
to supply the vacuum. The product is then sent to an oven where it is dried and sent to the manufacturer. It is the chemical engineer’s task to
create a consistent incidence graph in which the valves in the diagram do not depend on each other’s feedback. An example of this is below:

V1, TI, PI, TI (+,+,+,+)
V2, TI, PI, TI (+,-,-,-)
V3, PI, TI, PI (+,-,-,-)
Based on the table above, the incidence graph is consistent because the valves are not dependent on one another for feedback.
2.3.3.2: INCONSISTENT GRAPHS

Inconsistent pathways happen when pathways from the same starting node, that lead to the same target node, give different output.
For example,
Pathways from node 1 to target node 4 do not yield the same output. These are called inconsistent pathways.

1,4 (+,+)
1,2,4 (+,+,-)
1,3,4 (+,-,+)
From the table and graph above, all pathways from node 1 to node 4 do not yield the same outcome. Node 4 is being decreased while
simultaneously being increased.
An inconsistent graph is made up of only inconsistent pathways for all possible pathways between any combinations of two nodes. The
graph below is an example of an inconsistent graph.

1,2 (+,+)
1,3,2 (+,-,-)
1,2,1 (+,+,-)
1,3,2,1 (+,-,-,+)

1,2 : (+,+)
1,3,2 : (+,-,-)
Since the two pathways cause different changes to node 2, this sub-pathway is inconsistent.
1,2,1 : (+,+,-)
1,3,2,1: (+,-,-,+)
Because all sub-pathways in this incidence graph are inconsistent, the incidence graph is inconsistent.
2.3.3.3: PARTIALLY CONSISTENT GRAPHS

Partially consistent graphs are made up of both consistent and inconsistent pathways. If you increase node 1, all paths leading to node 4
cause a decrease in the node which means it is a consistent pathway. If node 1 is increased, one path leading to node 2 causes an increase
and another path causes a decrease in the node. This is an inconsistent pathway. Since this model contains both at least one consistent and
one inconsistent pathway, it is a partially consistent model.

1,2,3,4 (+,+,+,-)
1,2,4 (+,+,-)
1,4 (+,-)
1,2 (+,+)
1,4,2 (+,-,-)

1,2,3,4: (+,+,+,-)
1,2,4 : (+,+,-)
1,4 : (+,-)
Since the three pathways cause the same change to node 4, this sub-pathway is consistent.
All pathways leading from node 1 and to node 2 -
1,2 : (+,+)
1,4,2 : (+,-,-)
In this graph there are sub-pathways that are consistent and inconsistent, therefore the incidence graph is partially consistent.
2.3.3.4: SUMMARY
Just because a process is inconsistent does not mean that the desired results cannot be achieved. All it means is that you cannot consistently
achieve the desired results. Take the example of the flow of a final product. In a consistent process, no matter which path is taken, the flow
of the final product will always respond the same, whether desirable or not. In an inconsistent process, changing the path will not always
cause the flow of the final product to change the same way, but rather depends on the path taken.
In all, inconsistent processes are not bad. Often, processes will be very complex, containing many variables, and inconsistent processes will
be unavoidable. Care must just be taken with complex, inconsistent processes to achieve desirable results.
You are a process engineer and you have been assigned to devise a process control for a new reactor. The optimum operating condition for
the reaction inside the reactor is 220 oC ± 10oC. The reaction is exothermic and often goes way past the optimum temperature. It is
imperative that the reaction should not go above 250oC. The reactor has a built in heater (H1) to keep temperature at its optimum
temperature and a cooling water feed to cool the reactor in case temperature goes past 250oC. The reactor also has temperature sensor
installed to monitor the temperature. If temperature goes above 250oC, the desired mechanism is to increase the flow of water in the cooling
jacket surrounding the reactor and to turn off the heating element inside the reactor until it reaches its optimum temperature. As a control
redundancy mechanism, a pressure sensor is also installed inside the reactor. The pressure control will send signal to the cold water feed if
pressure reaches 1.2 atm which is the critical pressure. Draw a causation graph governing the situation when the temperature gets too hot
inside the reactor. Base your graph on the following variables:
Note: (There may be more than one correct solution to the assigned problem)
a. temperature
b. temperature sensor (T1)
c. cooling water feed (V3)
d. heating element (H1)
e. pressure
f. pressure sensor (P1)
High temperature will 'activate' the temperature sensor, telling it that the temperature has exceeded the optimal temperature and an action
need to be done to bring it back to the desired temperature. The temperature sensor would then open the valve for the cooling water feed
which would cause the water to flow through the jacket. This results in the temperature inside the reactor to decrease. The reverse happens
for the heating element. An activation of the temperature sensor would translate to an inhibition of the heating element in which the heating
element would stop increasing the temperature of the reactor. Temperature and pressure are directly related by ideal gas law and thus affects
each other in both ways. Thus, for the control redundancy mechanism, an increase in pressure would activate the pressure sensor which
would then activate the cooling water feed to cool the temperature for the reactor.
Logic control:
IF T1 =< 220 oC, THEN H1 is on, ELSE H1 is off
IF T1 >= 250 oC, THEN V2 and V3 is open, ELSE V2 and V3 is close
IF P1 >= 1.2 atm, THEN V2 and V3 is open, ELSE V2 and V3 is close

After completing the assignment of devising the temperature control, your boss asked you to devise a mechanism to control the flow of raw
materials and products. Draw a causation graph showing negative feedback relations between L1, V1 and V4 when the level inside the
reactor is increasing (this is an example of an inconsistent pathway). The normal operating conditions are shown below:
V1 is 50 gal/min, V1 is usually open
V4 is 50 gal/min, V4 is usually open
L1 is 20 meters, L1 should not go above 30 meters and should not be below 10 meters
The level sensor monitors the fluid height inside the jacketed CSTR. Assuming the level inside the reactor is increasing steadily, the level
sensor will send a negative feedback signal to V1 to decrease the input flow.
Logic control:
IF L1 > 30 meters, THEN V1 is closed
IF L1 < 10 meters, THEN V4 is closed

Determine whether the following causation graph is consistent, inconsistent or partially consistent. Include your reasoning.
Solution:
All relationships are consistent. Each variable 1,2,3,4 has a simple relationship which is monotonic among the other variables in the set.
{1,2} (+,-)
{1,2,3} (+,-,-)
{1,3} (+,-)
{1,4,3} (+,+,-)
{1,4} (+,+)
{2,1} (+,-)
{2,3} (+,+)
{2,1,3} (+,-,+)
{2,1,4} (+,-,-)
{4,3} (+,-)
2.3.7: REFERENCES
Woolf, Peter (2006). A Mini-Guide to Consistent Graphs.
Woolf, Peter (2006). Terms Review.
Angeli, David & Sontag, Eduardo D. (2003). Monotone Control Systems.
2.3.8: CONTRIBUTORS AND ATTRIBUTIONS

Authors: (September 26, 2006) Christopher Garcia, Anwar Stephens, Winardi Kusumaatmaja, Meng Yang Ng
Stewards: (September 5, 2007) Alexander Voice, Andrew Wilkins, Ibrahim Oraiqat, Rohan Parambi
This page titled 2.3: Incidence Graphs -Interpretations, Consistency, and Inconsistency is shared under a CC BY 3.0 license and was authored, remixed,
and/or curated by Peter Woolf et al. via source content that was edited to the style and standards of the LibreTexts platform; a detailed edit history is
available upon request.
2.4: EXCEL MODELING - LOGICAL MODELS, OPTIMIZATION WITH SOLVER
FOR NONLINEAR REGRESSION, SAMPLING RANDOM NUMBERS
2.4.1: INTRODUCTION
Microsoft Excel program is one of the most popular and useful computer programs for a wide variety of numerical applications. Excel
features several different functions, interfaces and graphing tools which can be applied in many fields. The following functions might be
especially useful for logical programming in Excel:
Logical functions (ex: IF, OR and AND): These functions can be used in control analysis, particularly in safety regulation (for example,
if the temperature exceeds X degrees, shut down the reactor and send cooling water to the jacket)
Solver: This function can be used to maximize, minimize or try to obtain an input value in a cell by varying referenced cells.
Random Number Generator. This can be used as a tool for probability or randomly selecting data points for analysis from a large pool of
points.
While it might be less accurate than other modeling programs, the Excel model provides good estimates, and facilitates observation and
understanding of process behavior while requiring minimum computer knowledge. The following sections will provide you with a better
understanding of these common Excel functions, followed by examples in engineering applications.
2.4.1.1: LOGICAL METHODS

A logical test is any expression that can be evaluated as TRUE or FALSE with specified input conditions. Commonly used logical tests are
IF, OR and AND statements. All three statements can be used in conjunction with each other for more complex and complete modeling.
Note: all functions and formulas in Excel require an equal (=) sign in front of the syntax.
2.4.1.1.1: EXCEL'S IF STATEMENT
The IF function in Excel tests values and formulas in the target cells and returns one specified value "if" the condition input evaluates to be
TRUE and another value if it evaluates to be FALSE. The IF statement in Excel is called out with the following syntax:
IF(logical_test,value_if_true,value_if_false)
Logical_test: any expression or correlation of values that can be evaluated as TRUE or FALSE.
Value_if_true: the value or expression that is returned IF the condition in the logical test returns TRUE.
Value_if_false: the value or expression that is returned IF the condition in the logical test returns FALSE.
The value_if_true and the value_if_false part in the syntax may vary and can be numerical values, expressions, models, and formulas. If
either of the values is left blank, the test will return 0 (zero) correspondingly.
Sample coding:
IF(1+1=2, "Correct!", "Incorrect!")
For the logical test part, 1+1 is equal to 2, and the iteration will return TRUE. The value_if_true is called for the TRUE response, and
"Correct!" will be output in the corresponding cell.
It is possible to nest multiple IF functions (called Nested IF) within one Excel formula. Such a formula is useful in returning multiple
answers, compared to a a single IF function which can return only two possible answers.
A Nested IF statement has the following syntax:
IF( logical_test1, value_if_true1, IF(logical_test2, value_if_true2, value_if_false) )
The following example shows how a Nested IF function is used to make a spreadsheet that assigns a letter grade to different scores.
Suppose that grade distribution is:
@ A= Greater than or equal to 80
@ B= Greater than or equal to 60, but less than 80
@ C= Less than 60
This logic function can be written in several different ways. If the score is tested to be the letter grade C first, then checked to see if it is a B
or an A, the syntax could be written as follows:
IF(A5<60,"C",IF(A5<80,"B","A"))
Note that the Nested IF statement can be placed in either the value_if_false OR the value_if_true positions. In the example above, it was
placed in the value_if_false position because of how the conditional was written.
Additional letter grades could be added by including more nested IF functions inside the existing function. One formula can include up to 7
nested IFs, which means you can return up to 8 different results.
2.4.1.1.2: EXCEL'S OR STATEMENT
The OR statement in Excel is similar to the IF statement. It returns TRUE if any argument or expression of values is TRUE, and it returns
FALSE if all arguments or expression of values are FALSE. The syntax for OR is:
OR(logical1,logical2,...)
logical1: there can be up to 30 arguments for the logical expression to test and return either TRUE or FALSE. If all of the arguments are
false, the statement returns FALSE for the entire expression.
Sample coding:
OR(1+1=2, 2+2=5)
In the logical1 part, 1+1 is equal to 2 and will return TRUE. In the logical2 part, 2+2 is not equal to 5 and will return FALSE. However,
since one of the test returns TRUE (logical1), the function returns TRUE even if logical2 is FALSE. TRUE will be output in the
corresponding function cell.
2.4.1.1.3: EXCEL'S AND STATEMENT
The AND statement returns TRUE if all its arguments or values are TRUE and returns FALSE if one or more argument or value is FALSE.
The AND statement syntax is:
AND(logical1,logical2,...)
logical1: similar to OR statement with up to 30 arguments to test that could return either TRUE or FALSE.
Sample coding:
AND(1+1=2, 2+2=5)
Logical1 will return TRUE for the correct calculation while logical2 will return FALSE for the incorrect calculation. However, since one of
the expressions is FALSE, the AND logical test will return FALSE even if one of the tests is true.
It is common and useful to understand how to implement different logical functions into each other. A good example would be:
IF(AND(2+2=4,3+3=5),"Right Calculation","Wrong Calculation")
During iteration, the AND statement should return FALSE since one of the expression/calculation is incorrect. The FALSE return will
output "Wrong Calculation" for the IF statement iteration.
2.4.1.1.4: OTHER USEFUL EXCEL FUNCTIONS
Excel has many useful functions that can help any engineer who is using Excel to model or control a system. An engineer may want to
know the average flow rate, number of times the liquid in a tank reaches a certain level, or maybe the number of time steps that a valve is
completely open or completely closed. All of these situations and many more can be addressed by many of the simple built-in functions of
Excel. Below are some of these functions and a brief definition of what they compute and how to use them. They are followed by an image
of a sample excel file with an example of each function in use. To see a list of all the built-in functions Excel has to offer, click shift-enter
on the cell.
AVERAGE(number1,number2,...) - Returns the average or arithmetic mean of the arguments. The arguments can be numbers, arrays, or
references that contain numbers.
CEILING(number1, significance) - Rounds a number up to the nearest integer or multiple of significance.
CHIDIST(x, degrees_freedom) - Returns the probability of the chi-squared distribution.
COUNT(value1,value2,...) - Counts the number of cells that contains numbers.
COUNTIF(range, criteria) - Counts the number of cells within the range that meet the given criteria.
COUNTIFS(criteria_range1,criteria1,criteria_range2,criteria2,...) - Counts the number of cells within multiple ranges that meet multiple
criteria.
FLOOR(number1,significance) - Rounds a number down toward zero, to the nearest multiple of significance.
FREQUENCY(data_array,bins_array) - Calculates how often a value occurs within a range of values (bins). The function will return a
vertical array of numbers having one more element than the bins_array. Once the formula is entered into the target cell, highlight the target
cell making a vertical array one cell larger than bins_array. Once selected press F2 to open the formula of the target cell and then press
CTRL+SHIFT+ENTER to execute your function.
LARGE(array,k) - Returns the k-th largest value in a data set. For example, the 2nd largest number.
MAX(number1,number2,...) - Returns the largest number in a set of values.
MEDIAN(number1,number2,...) - Returns the median of a set of values.
MIN(number1,number2,...) - Returns the minimum value in a set of values.
MODE(number1,number2,...) - Returns the most frequently occurring value in an array or range of values.
PERCENTILE(array,k) - Returns the k-th percentile of values in a range
ROUND(number1, num_digits) - Rounds a number to the specified number of digits.
SMALL(array,k) - Returns the k-th smallest value in a data set. For example, the 7th smallest number.
STDEV(number1,number2,...) - Returns the standard deviation of a given set of values ignoring text and logical values.
SUM(number1,number2,...) - Adds all of the values in a range of cells.
SUMIF(range,criteria,range,...) - Adds values of cells in a range meeting criteria.
VAR(number1,number2,...) - Estimates the variance of a given set of values.
2.4.1.1.5: NESTING
Though not complex or difficult, placing logical statements within one another can be a powerful practice. This is referred to as nesting.
An example demonstrating the use of a simple IF statement nested within another IF statement is shown below:
If one were to ask three individuals whether or not they agree that energy has a velocity, these individual will hypothetically respond in one
of three ways: "Yes","No", or (other). In this case, if the individual were to respond "Yes", then the questioner would respond "I agree!". If
the individual were to respond "No", the questioner would bitterly respond with a "Bah". And if the individual were to respond with a
nonsensical answer, then the questioner would respond "what .." in confusion.
The important point to get out of this example is that nested statements gives us the potential to process more information in one cell (or
line of programming code).
2.4.1.2: EXCEL'S DATA ANALYSIS TOOLS
Excel has several useful built in add-ins that can make data analysis much easier. Of these we will cover the solver tool, as well as the
ANOVA analysis.
2.4.1.2.1: SOLVER TOOL
The solver in Excel is a built-in tool with multiple capabilities. It could be used to optimize (find the maximum or the minimum value) a
specific cell (the target cell) by varying the values in other cells (the adjustable cells). It may also be used to solve a system of non-linear
equations. Solver formulas embedded in the spreadsheet tie the value in the target cell to the values in the adjustable cells. The user can also
specify numerical range constraints for the cells involved.
Optimization Model
In optimization, one seeks to maximize or minimize the value of a real function. Just as a student in a calculus class might use optimization
techniques to find a local maxima or minima of a complex function, a CEO might use similar techniques to find how much product their
company should manufacture to maximize profits. There are numerous optimization techniques, many specialized for specific types
problems. For the most part, MS Excel uses gradient-based techniques.
The solver tool is located in the Tools menu. If it does not appear there, it must be added in by selecting Add-ins and selecting the
appropriate check box.
The solver window contains the following fields:
Set Target Cell: This is the cell containing the value you want to optimize.
Equal To: Choose whether to maximize the value, minimize the value, or set it to a specific number.
By Changing Cells: Here you specify all the cells that can be varied when optimizing the target cell value.
Subject to the Constraints: Constraints are optional, but can be added by clicking Add and defining each constraint.
List of Available Constraints
<= (Less than or equal to) Stipulates that a selected cell must be less than or equal to a certain value.
= (Equal to) Stipulates that a selected cell must be equal to a certain value.
>= (Greater than or equal to) Stipulates that a selected cell must be greater than or equal to a certain value.
int (Integer) Stipulates that a selected cell must be an integer
bin (Binary) Stipulates that a selected cell must be equal to 1 or 0.
The window where these constraints can be entered is shown below.
It is important to make sure that constraints do not violate each other otherwise solver will not work without reporting an error. For
example, if a cell must be great than 2 and less than 1 Solver will not change anything.
Once all the fields in the Solver window are completed appropriately, click Solve. If solver can find a solution a window will pop up telling
you so, and the values will appear in the target and adjustable cells. If solver cannot find a solution it will still display the values it found but
will state that it could not find a feasible solution. Be aware that Excel solver is not a powerful enough tool to find a solution to nonlinear
equations, where programs such as Matlab might be a better choice.
* Guided Example: Suppose you have an Excel spreadsheet set up as shown below, with cell C6 depending on cell B6 according to the
quadratic relationship shown, and you want to minimize the value in C6 by varying the value in B6.
1) Open up the solver window.
2) Input C6 in the set target cell field, check the min button, and input B6 in the by changing cells field:
3) Click solve to find the solution.

You can place further constraints or further specify how you want your task solved under "Options." Some of the most commonly used
specifications include:
Max Time: This is the maximum amount of time excel will spend trying to find a converging value before giving up.
Iterations: This is the number of iterations excel will perform to converge to a value.
Precision: This is related to how close excel has to get before a value is "acceptable." In general, making this value several orders of
magnitude smaller will yield a more accurate result, but will take more time to do so, and may have a harder time converging.
Tolerance: This is also related to how accurate a solution is, this is related to the percent error. Making this smaller should also yield a
more accurate result, but with more difficulty finding a solution.
Convergence: Controls the amount of relative change that occurs during the last five iterations that Solver performs. The smaller the
number of the convergence the less change that will occur.
Box Options Click on the box to select the following options:
Assume Linear Model: This will speed up the time to find a solution because it will make Solver assume that the model is linear.
Assume Non-Negative: Solver will make the lower limit of adjust zero for all the cells that are programmed to be adjusted in the model
that are not limited by a chosen constraint.
Use Automatic Scaling: Used for large differences in magnitude to scale the inputs and outputs.
Show Iteration Results: Solver will pause between each iteration to show the results of that iteration in the model.
Estimates Specifies to Solver on how to use initial estimates for iterations
Tangent: Uses linear extrapolation by utilizing a tanget vector.
Quadratic: Uses quadratic extrapolation. Best to use for nonlinear models.
Derivatives Specifies what type of differencing is used for partial derivatives for both the objective and constraints.
Forward: Constrained values change slowly between iterations.
2.4.1.3: NONLINEAR REGRESSION

Excel's solver function can also be used to find a solution for two-variable non-linear regression. The description of the data by a function is
carried out by the process of iterative nonlinear regression. This process minimizes the value of the squared sum of the difference between
the actual data and the predicted value (residual).
There are three basic requirements for using solver including:
1. Raw data
2. A model equation to find predicted values
3. Initial guesses for varying values
Before using the solver function, the raw data should be organized into separate columns for the independent variable (first column) and
dependent variable (second column). A third column will show the predicted value that is obtained by using the model equation mentioned
earlier, and must reference various single cells that have initial values that Excel will use as a starting point for the regression. A fourth
column should be created with the square of the difference between the predicted value and the actual experimental value (raw data); this is
also known as the square of the residual. A cell at the bottom of this fourth column should include the sum of all the squared difference
values above it. Two more cells anywhere else in the worksheet should be designated as the initial guesses for the values to be varied to find
the solution.
It is important to choose a reasonable value for the initial guess. This will not increase the chance that Excel arrives at a global optimization
rather than a local optimization, but also reduces the time it takes for Excel to solve the system of non-linear equations.
Now, solver can be applied to the spreadsheet to optimize the regression. Refer to section 3.1 for how to access and use solver. The bottom
most cell of the fourth column that has the sum of the squares of the residuals will be the target cell. Since the residuals are essentially the
error, the sum of the residuals squared should be as small as possible, so the min circle should be selected. The initial guesses should be
entered into the "by changing cells". Constraints can be entered if the values of the cells are limited to a certain range. Solver will produce
the values of the varying parameters that yield the best fit to the nonlinear data.
Guided Example: A screen shot of a simple population problem is shown below. The yellow highlighted cells are the initial guesses for
the variables. In this case, the initial population (P0, cell E2) is entered as a constraint.
2.4.1.3.1: ANOVA ANALYSIS

ANOVA is another tool on Excel that can be used for data analysis purposes. It is used primarily to compare continuous measurements to
determine if they are sampled from the same or different distributions. This is usually done on set of data that consists of numerous samples
and/or groups. Single factor ANOVA analyzes the difference between "groups" which can be different units, sensors, etc. On the other hand
2 way ANOVA analyzes the differences on the basis of both samples (trials) and groups, allowing you to compare, which one most affects
the system. Two-way ANOVA can be employed for two different cases: data with or data without replicates. Data without replicates is used
when collecting a single data point for a specified condition, while data with replicates is used when collecting multiple data points for a
specific condition.
The telling result of ANOVA analysis is a p-value telling you if there is a statistically significant difference between data sets. A more in
depth discussion of the uses and implementation of ANOVA can be seen here.[1]
2.4.1.4: RANDOM NUMBER SAMPLING

A useful numerical technique is the random number generator. This generator uses random number sampling to create a predictive model.
This can be especially useful in probability problems.
A typical random number generator produces any number between 0 and 1 with equal frequency. At each step a random number is
generated that picks a possible outcome. A typical Excel statement would appear as the following: IF (RAND()<X, value1, value2). This
means that X percent of the time, you would get a value1 and 1-X percent of the time, you would get a value2. For example, suppose a
situation exists with two possible outcomes, A and B, where A takes place 70% of the time. If a random number is generated that is greater
than 0.3, then the outcome is A. Inversely if the random number is less than or equal to 0.3, then B is the outcome. This would be
represented as IF(RAND()<0.7, A, B). Each step within the generator does not affect the next step, i.e. each step is determined individually.
A process in which previous outcomes do not affect later outcomes is known as a Stochastic method or "Monte Carlo" method. A Monte
Carlo simulation uses a model that takes random input that follows an assumed distribution (in this case, 70% A and 30% B) and produces
output that models some phenomenon. The following is a simple example of such a simulation:
The Excel is equipped with a random number generator which can be called out using the RAND function. A random number between 0
and 1 can be generated in a cell by the following command:
In Microsoft Excel 2007 the RAND function returns a new random number every time the worksheet is calculated or when F9 is pressed. If
the same random number is needed to be referred to more than once in an Excel IF, AND, OR statement, one can create a column of random
numbers using =RAND(), and refer back to these cells in an adjacent column.
There are many ways to manipulate the random numbers that Excel produces. A random number can be generated within a different range
by multiplying the RAND function by a constant, i.e. multiply by 100 to change a decimal to percent. Another option to generate numbers
within a different range is the RANDBETWEEN function. RANDBETWEEN generates a random integer within a specified range. The
following command generates a random integer between -6 and 21:
It is also possible to modify your random numbers by nesting the RAND function within another function. For example, one can square the
results of one distribution to create a different distribution of numbers.
Another way to utilize the RAND function is to take a random sampling of a data set. If one has a spreadsheet with several hundred data
points, and would like to analyze only 100 random points, another column can be created titled "Random Number" or something similar. In
the first cell of this column, type the =RAND() function. Then drag this function all the way to the end of the data set. At this point, the user
can go to the DATA drag-down menu and select Sort... In this menu, select your random number column in the SORT BY drag down menu.
Either ascending or descending can be selected, as desired. Once the spreadsheet is sorted by the random numbered cells, the first 100 data
points can be used as a statistically random cross-section of the data.
The sort tool can be useful in cases when you are trying to find out the probability of having certain values above a certain threshold. For
example, if the 100 data points that you have represent the temperature variations occurring within a system, you can use sort tool to
determine how many temperature values exceed your threshold. From that, you can calculate approximately the probability of the system
temperature exceeding the threshold. Two screenshots of the sort process are shown below. Microsoft Excel 2007 uses a different menu,
allowing for one to choose between sorting from Largest to Smallest or vice versa. A screenshot below shows this new menu.
Sorting Setup:
Results to Randomization:
Microsoft Excel 2007:
2.4.1.5: WORKED OUT EXAMPLE 1

The Arrhenius equation defines the relationship between a reaction rate constant k and temperature:
−E /RT
k(T ) = Ae
where T is the absolute temperature, A is the frequency factor, E is the activation energy, and R is the universal gas constant.
It is frequently desirable to be able to predict reaction rates at a given temperature, but first you need to know the values of A and E. One
way to obtain these values is to determine the reaction rate constant at a few temperatures and then perform a nonlinear regression to fit the
data to the Arrhenius equation. Given the following rate constant data for the decomposition of benzene diazonium chloride, determine the
frequency factor and activation energy for the reaction.
−1
(s ) (K)
0.00043 313.0
0.00103 319.0
0.00180 323.0
0.00355 328.0
0.00717 333.0
Solution: The following Excel file contains a solution to the problem that uses Excel's solver tool to perform the nonlinear regression and
determine the values of A and E:
Example 1
The spreadsheet is set up as described in the nonlinear regression section above.
1. Open solver (Tools / Solver)
2. Set the Sum of Residuals Squared value cell as the target cell
3. Click the Min radio button
4. Set Pre-exponential Factor and Activation Energy as the adjustable cells. For this problem, keep A between 1E+13 and 1E+14 s^-1 and E
between 9.5E+4 and 1.05E+5 J
5. Sometimes, depending on what version of excel you are using, a message will come up like the one below. If this happens, click OK. The
graph will update with the new curves.
6. Click Solve and observe the values of A and E calculated by Solver and the changes to the plot showing the rate constant as a function of
temperature.
NOTE: If solver is not fitting the data well, click the options button in the solver window and make the precision value smaller by several
orders of magnitude.
2.4.1.6: WORKED OUT EXAMPLE 2

Kinetics can be used to determine the probability of a reaction taking place between gas-phase molecules. Consider molecules A and B,
which have a 60% chance of reacting upon collision. Create an Excel spreadsheet that utilizes a random number generator and an if-
statement to model an individual collision.
Solution
In this Excel file, sample coding for Random number generator and if-statement is shown.
Collision_ex
Random number generator is used to output numbers between 0 and 1 for four trials, ten trials, and one hundred trials on sheets one, two,
and three, respectively. Formula for the random number generator is employed in the cell where you want the number to be output (in this
case, B6 is the cell). An IF logical test is then employed to test whether the number follows the requirement or not. The IF-statement cell
will then output the corresponding result based on the logic test. In this collision case, if the random number function generated a number
greater than 0.4, the logic test will return TRUE and then output "Reacts" as the result. Assuming that the random number generator
produces an even distribution between 0 and 1, the random number will be greater than 0.4 sixty percent of the time. Coding of the IF-
statement is shown:
On the other hand, if the number generated is less than 0.4, the logic test returns FALSE and "Does not React" will be shown as the result.
Different numbers are generated each time the spreadsheet is updated, thus different result will return for the IF function. You can also
modify the numbers and return statement to better understand the operation of the functions.
The number of reactions is summed up and divided by the total number of trials to compare the predicted reaction probability with the value
given in the problem. It can be seen by comparing sheets one, two, and three that increasing the number of trials decreases the variation in
the predicted reaction probability.
2.4.1.7: ACCESSING EXCEL TOOLS ON WINDOWS

On many computers, excel tools such as solver and the data analysis package are not already included and must be installed before used.
This section will go over how to access these add-ins for those using Windows opposed to Macs.
1) Click on the windows symbol in Excel 2007. At the bottom of the opened window, select the Excel options box.
2) On the left side panel you will see a section labeled "Add-ins". Select this.
3) Here you will see a list of "Inactive Applications" that are available for excel. Of these the ones we are interested in are the "Solver Add-
In" and the "Analysis ToolPak." Select one of these and click "Go" at the bottom of the window.
4) A new window will appear asking you which components of the package you would like to install. Select the needed tools (solver and
analysis toolpak). Click "Ok" and windows should install the corresponding tools into excel.
5) Once done, these tools will be found under the data section of Excel, on the far right of the window.
2.4.1.8: REFERENCES
1. =RAND()
=RANDBETWEEEN(-6,21)
Collision Example Excel File:
=IF(B6>=0.4,"Reacts","Does not React")
Bender, E.A. An Introduction to Mathematical Modeling, Mineola, NY: Dover Publications, Inc.
Fogler, H.S. (2006). Elements of Chemical Reaction Engineering, Upper Saddle River, NJ: Prentice Hall Professional Technical
Reference. ISBN 0-13-047394-4
Microsoft, Excel, Proprietary EULA, www.microsoft.com.

Authors: (12 September 2006 / Date Revised: 19 September 2006) Curt Longcore, Ben Van Kuiken, Jeffrey Carey, Angela Yeung
Stewards: (13 September 2007) So Hyun Ahn, Kyle Goszyk, Michael Peterson, Samuel Seo
This page titled 2.4: Excel Modeling - logical models, optimization with solver for nonlinear regression, sampling random numbers is shared under a CC
BY 3.0 license and was authored, remixed, and/or curated by Peter Woolf et al. via source content that was edited to the style and standards of the
LibreTexts platform; a detailed edit history is available upon request.
2.5: NOISE MODELING - WHITE, PINK, AND BROWN NOISE, POPS AND
CRACKLES
2.5.1: INTRODUCTION
Noise is all around us in all sorts of forms. In the common use of the word, noise refers to sound. However, noise can more accurately be
thought of as random variations that are always present in one or more parts of any entity such as voltage, current, or even data. Rather than
thinking of noise only as an acoustic term, it should be thought of more as a random signal. Noise can be the inherent fluctuations in some
part in a system (ie. temperature at a given point) or it can be the unavoidable interference on a measurement from outside sources (ie.
vibrations from a nearby generator blur measurements from a pressure transducer). The static interference on your radio, the ‘snow’ on your
television, and the unresolved peaks on an infrared spectroscopy report are all examples of noise.
Chemical engineers can use statistical properties to characterize noise so they can understand a current process or develop an optimal
process. By characterizing the noise and determining its source, the engineer can devise methods to account for the noise, control a process,
or predict the path of a system. For example, when chemical engineers design plants, they use characterizations of noise to determine the
best control scheme for each process. Mathematical modeling can be used to characterize and predict the noise of a given system. In
modeling, to simplify the representation of general trends that reoccur, noise is classified in two major categories: frequency based and non-
frequency based. Frequency based noise consists of the colors of noise, and non-frequency based noise includes pops, snaps and crackle.
The purpose of the following sections is to give you a qualitative understanding of the sources of noise and how it can be characterized and
handled. Examples will be given to help you gain a quantitative understanding of how noise relates to controlling and modeling chemical
engineering processes.
2.5.2: GENERAL TYPES OF NOISE

Three general types of noise can be categorized in describing a process, and are as follows:
1) Gaussian noise is usually not dependent on time, meaning that it is random and not systematically planned. The amplitude of the
frequency can vary, making a crackling notation or sound (see #Crackles below). Some examples of Gaussian noise can include splashing in
a tank or unplanned interruption in a sensing device.
2) Drift noise is correlated to time, and has random movement. Examples of drift noise can include fouling or catalyst decay, in which the
potency of the substance may decline over time as decaying occurs. Stock price fluctuations also model somewhat like drift noise.
3) Shot noise may be defined as sporadic and short bursts of noise, in which the amplitude is similar among the bursts of noise. Shot noise
can be correlated to pops (see #Pops below), in which at random times, the same shot noise is witnessed with the same amplitude. Examples
of shot noise include partial clogging or jamming in the process in which the same amplitude will be seen by the noise whenever the
clogging or jamming in the process occurs. Another example is customer demand of the product. If the control system (when trying to
optimize it) depends on the customer order, and customer demand is not consistent at all times (meaning downtimes for orders) but the order
amount is the same when orders are placed, then it will be effected by the shot noise described by the customer order (customer demand).
Microsoft Excel Formulas for 3 Types of Noise

Noise Type Excel Formula
Gaussian Noise = [Original value] + RAND() - RAND() + RAND() - RAND()1
Drift Noise = [Previous value] + RAND() - RAND() + RAND() - RAND()1
Shot Noise = if(RAND()>0.9, [Original value]+[single error], [Original value])
1
The more add/subtract cycles of RAND() you use, the closer approximation you'll get for a Gaussian distribution
of random errors.
As you would see these in an Excel file:
\ A B C D
1 Time Gaussian Noise Drift Noise Shot Noise
2 1 original value original value original value
3 2 = $B$2 + RAND() - RAND() + RAND() - RAND() = C2 + RAND() - RAND() + RAND() - RAND() = if(RAND()>0.9, $D$2+[single error value], $D$2)
4 3 = $B$2 + RAND() - RAND() + RAND() - RAND() = C3 + RAND() - RAND() + RAND() - RAND() = if(RAND()>0.9, $D$2+[single error value], $D$2)
Note: These formulas can be extended for any length of time. Longer periods of time are useful for visualizing drift noise.
2.5.2.1: COMBINED TYPES OF NOISE
Sensors, processes, demands, etc. often do not behave with simple noise models. However, most noise can be derived from the three general
types of noise. An example of a process with random and shot noise is shown in the thumbnail.
This is a graph of the total daily electricity usage for a plant during a normal operation period. As we can see, there are minor fluctuations
occurring every day, but some days there are large changes of electricity usage. These large changes can be seen as shot noise when the
consumption is viewed on this daily scale.
Possible reasons for the minor fluctuations could be due to electric heaters kicking on and off. Changes in the operation of smaller pumps
can cause these changes in electricity demand.
The large shots could be due to an energy-intensive part of the process that only operates when needed. The positive shots could be due to
this process kicking on and the negative shots could be due to a large process kicking off when not needed.
With this graph, we can say that there is shot noise present with the large changes in consumption, but there is also random noise present
with the minor daily fluctuations. We would not say there is drift noise present because the consumption remains around a constant mean
value and does not drift away from there.
However, the total consumption when viewed from a different scale may show other forms of noise. If we viewed the hourly usage of
electricity on one of the days with shot noise, it may resemble drift noise on the hourly scale.
2.5.3: COLORS OF NOISE

The first category of noise is frequency-based noise, classified by the colors of noise. In understanding the colors of noise, we need to
understand the process of converting a given set of data to a point where it can be classified as some color of noise. The section below,
“Modeling Noise,” goes into more detail as to how this is done. At this point, it is important to understand the main colors of noise which
are possible in a system. Noise is classified by the spectral density, which is proportional to the reciprocal of frequency (f ) raised to the
power of beta. The power spectral density (watts per hertz) illustrates how the power (watts) or intensity (watts per square meter) of a signal
varies with frequency (hertz).
1
P SD ∝
β
f
where β ≥ 0.
Just like colors of light are discriminated by their frequencies, there are different colors of noise. Each beta value corresponds to a different
color of noise. For example, when beta equals zero, noise is described as white and when beta equals two, noise is described as brown.
2.5.3.1: WHITE NOISE

One major type of noise is white noise. A defining characteristic of this type of noise is it has a flat power spectral density, meaning it has
equal power at any frequency. β equals zero for white noise.
White noise and white light have similar properties. White light consists of all visible colors in the spectrum while white noise is created by
combining sounds at all different frequencies. Pure white noise across all frequencies cannot physically exist. It requires an infinite amount
of energy and all known energy is finite. White noise can only be created within a specific and defined range of frequencies. Again making
an analogy to white light, for a small band of frequencies, visible white light has a flat frequency spectrum. Another property of white noise
is that it is independent of time. This factor also contributes to the idea that white noise cannot physically exist since it is impossible to have
noise be completely independent of time.
A sample power spectral density chart of white noise is shown below. As can be seen from the chart, the power over all displayed
frequencies is essentially the same, suggesting the PSD is equal to 1.
image from en.Wikipedia.org/wiki/White_noise on 9/18/2006

It is important to note that white noise is not always Gaussian noise. Gaussian noise means the probability density function of the noise has
a Gaussian distribution, which basically defines the probability of the signal having a certain value. Whereas white noise simply means that
the signal power is distributed equally over time. For Gaussian white noise one finds the number of measurements greater than the power
spectral line as measurements that are less than the power spectral line with majority falling at the spectral line. Gaussian white noise is
often used as a model for background noise in satellite communication. White noise can also come from other distributions, such at the
Poisson distribution. Poissonian white noise will look like a normal distribution that has been shifted to the left for small number of
measurements. For larger number of measurements Poissonian white noise will look like the normal distribution (also known as gaussian
white noise). A Poissonian white noise model is useful for systems where things happen in discrete amounts such as photons being
transmitted or the number of cars that go through an intersection.
2.5.3.2: PINK NOISE

Pink noise is a signal whose power spectral density decreases proportionally to the inverse of the frequency, where the β value is equal to
one. Because of its decrease in power at lower frequencies when compared to white noise, pink noise often sounds softer and damper than
white noise.
Pink noise is actually what is commonly considered audible white noise and is evident in electrical circuits. Pink noise becomes important
to measure in things such as carbon composition resistors, where it is found as the excess noise above the produced thermal noise. Below is
a power spectral density chart for Pink noise.
image from en.Wikipedia.org/wiki/Pink_noise on 9/18/2006
2.5.3.3: BROWN NOISE

When β equals 2, the noise is Brownian. Brown noise is signal noise created by Brownian, or random, motion. It is a form of unavoidable
interference in the transmission of information and can be compared to a random walk that does not have a clearly patterned path. Brown
noise has more energy at lower frequencies. Compared to white or pink noise, brown noise is even more soft and damped. It can be
generated by integrating white noise. Interference from outside sources, such as vibrations from nearby machinery or background light, in
instrument readings usually have a brown noise pattern.
Below is a power spectral density chart for brown noise. From the charts of brown noise and pink noise, it can be observed that brown noise
loses power as frequency increases at a much faster rate than that of pink noise.
image from en.Wikipedia.org/wiki/Brown_noise on 9/18/2006
2.5.3.4: BLUE NOISE

Blue noise is a signal whose power spectral density increases proportionally to the frequency, where the β value is equal to negative one. An
example of would be a reversible, exothermic, batch reaction where the products are being accumulated at an overall constant rate, but the
reverse reaction occurs when the temperature gets too high. Thus it is overall increasing in product, but has small fluctuations as it
increases. Below is a power spectral density chart for blue noise.
image from en.Wikipedia.org/wiki/Blue_noise#Blue_noise on 12/11/2008
2.5.3.5: PURPLE NOISE

When β equals -2, the noise is violet. Purple noise has more energy at higher frequencies. It can be generated by differentiating white noise.
Below is a power spectral density chart for purple noise. From the charts of purple noise and blue noise, it can be observed that purple noise
gains power as frequency increases at a much faster rate than that of blue noise.
2.5.3.6: MODELING COLORS OF NOISE

Characterizing noisy signals is important to a chemical engineer so that he or she can determine the sources of noise and how to account for
it. Sometimes, the noise will be a characteristic of the system. In these cases, the system can not have a very sensitive control because if the
sensor for the controlling device is too sensitive, it will respond to the noise of your system and end up controlling your system in a noisy
way. Other times, the noise will be from unavoidable disturbances to the system. In these cases, the noise can be damped or accounted for in
an appropriate manner.
In order to apply the knowledge of the colors of noise we need to understand the process of converting a given set of data to a point where it
can be classified as some color of noise. The general process of classifying noise follows.
Data Curve to Fit Data Fourier Transform Power Spectral Density Classification of Noise Color
The classification of noise color becomes important when we need to have an estimate of future data. In other words, given the trends in our
current data we can use it to estimate future data. This is done as follows,
Noise Color Power Spectral Density Fourier Transform Curve Generate Data From Curve
The power of a noise signal is detected at a certain frequency. Then a plot of the log(power) vs. the log(frequency) can be constructed, and
the slope of the line gives the beta value. Following a backward thought process, one can produce a certain color of noise by creating
frequency components which have a value generated by a Gaussian distribution and then scaling by the appropriate beta power of
frequency.
A general method to characterize and model noise is explained below.
1. DATA
Data is what the signal transmits. The signal is dependent on what you are measuring or modeling. The data can be collected, for example, if
you're measuring temperature in a reactor, then your data is the temperature readings from a thermocouple at a certain position in the reactor
over a period of time.
2. CURVES
After you collect the data, you plot the data and find a best fit equation, , for that set of data. A math program can be used to find a best
fit equation. Microsoft Excel can be used for a simple model and programs such as Polymath or Mathematica can be used for more complex
models. A simple equation could be in the form of . The coefficients A and B can be varied to fit the data. Additional terms can
be added for more accurate modeling of complex data. The equation you find can then be used to predict and model future signal noise.
3. FOURIER TRANSFORMS
A Fourier transform is a tool used to convert your data to a function of . In this form, the noise can be more easily characterized. Apply a
Fourier transform to the curve, , you fit to your data to generate a relation for the power spectral density. The Fourier transform can be
performed by a computer depending on the complexity of (or see "Simplifying the Fourier Transform" below). The transform is the
integral shown below.
∞
−jwt
X(w) = ∫ x(t)e dt
−∞
where;
is the equation of the curve to fit the data.
is the exponential form of writing the relation . The j in the second term here indicates that that term is
imaginary.
is the frequency
4. POWER SPECTRAL DENSITY
This value is attained by simplifying and squaring the Fourier Transform. Since the series of simplifications needed are fairly complex, the
resulting power spectral density equation is displayed below.
−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−
2 2
P S D = √ (∑ x(t) cos(wt)) + (∑ x(t) sin(wt))
At this point we attain a numerical value for the PSD at the particular frequency, . These numerical PSD values can be plotted versus
frequency to obtain the PSD chart.
5. CLASSIFICATION OF NOISE
The summation is repeated over and over for different frequencies, . A plot of the PSD vs. frequency, is made with these values. Once
this is done a best fit line is applied to the data points which will give you the characterization relation . Based on this we can
then classify the noise as a color of noise.
image from en.Wikipedia.org/wiki/Violet_noise#Violet_noise on 12/11/2008
THE REVERSE PROCESS
Knowing how to convert data to a color of noise is only half the problem. What if we know what type of noise is possible and we need data
from it for a given process? Knowing the noise color means that we know the power spectral density relation to the frequency. From here
onwards we follow the reverse route as that taken to get to the noise color by using the inverse Fourier transform instead.
Inverse Fourier Transform
∞
1
jwt
x(t) = ∫ X(w)e dt
2π −∞
where;
is the equation of the curve to produce future the data.
is the exponential form of writing the relation . The j in the second term here indicates that that term is
imaginary.
is the frequency
We can either use it to produce a curve or, by making a similar simplification as we did for the Fourier transform, to generate data directly.
This reverse process should be trivial to someone who worked through the forward process.
2.5.4: POPS, SNAPS & CRACKLES

The second category of noise is non-frequency based noise. Three examples of this non-frequency noise are pops, snaps and crackles. The
study of these noise types is fairly new, and not a lot is known about how to deal with them when they arise in our research or studies.
However, one can successfully distinguish and classify these noise types by noticing certain properties which are typical of them.
2.5.4.1: POPS
At one end of the extreme of non-frequency noise is what is defined as pops. Pops are infrequent random spikes in noise of approximately
the same amplitude. An example of this would be a loose wire that is usually completing the circuit. However, it randomly disconnects for a
split second and then recompletes the circuit for normal operation. Chemical engineering processes may inherently have pops that are
unpredictable. Safety features can be added to a process, such as pressure relief valves, to handle the pops when they do occur. The
following is an image of what pops might look like in a system:
2.5.4.2: SNAPS
On the other end of the non-frequency noise spectrum there are snaps. Snaps are single independent events that occur only once. When a
pencil is bent and it breaks under the force is an example of snapping. Another example would be the snapping of a piece of chalk or
exploding a bomb. The “popping” of one’s knuckles is really snapping since it is an independent event, unless you snap all your knuckles all
at once or one after the other really fast. Just like pops, snaps are unpredictable, and safety features should be added to the system to handle
the snaps. The following is an example of a snap:
2.5.4.3: CRACKLES
In between popping and snapping there is crackling. A very common example for crackling is the sound heard coming from a burning piece
of wood. Like popping, there is a non-frequency or irregularity in which the crackles occur. In addition, there is also an irregularity of the
amplitude of the crackle. In the case of the fire, not only can you not predict when the crackle sound would be heard, but you cannot predict
how loud it will be either. Furthermore, there is a universality condition associated with crackling. Regardless of the scale, the similar
randomness in repetition and amplitude should be observed.
In dealing with this universality condition the concept of a critical exponent arises. For example, if we are looking at the same crackling
effect, S , over a larger period of time the two would have to be equal after scaling the larger one.
T
⟨S ⟩small (T ) = A⟨S ⟩large ( )
B
for some scaling factors and .

where ⟨S ⟩ is S average.
If the time scale is expanded by a small factor B = ( 1
(1−δ)
) then the rescaling of the size will also be small, say 1 + aδ .
⟨S ⟩(T ) = (1 + aδ) ∗ ⟨S ⟩(1 + δT )
For small δ we have

d⟨S ⟩
a⟨S ⟩ = T
dT
Solving gives
a
⟨S ⟩(T ) = So T (2.5.1)
The exponent a is called the critical exponent and is a universal prediction of a given theory. There are other concepts which are used as a
check for universality for times when the critical exponent cannot be used but the most common one is the critical exponent. The following
is an example of what a crackle might look like in a system:
 EXAMPLE 2.5.1: CRACKLE
A lead engineer monitoring the instantaneous flow rate of a coolant used to cool an exothermic reactor gathers data set ‘a’ over a 30hr
period. Later on that week he gathers data set ‘b’ for the same coolant flow rate. For both sets of data he determined that the way the
instantaneous coolant flow rate is reacting to the exothermic reactor is optimal to the reaction process. For his records he would like to
represent the noise in the data in a compact form. He wants you to first characterize the noise and then provide any other information,
which will generalize the trend in the data. Once you have completed these tasks, report back to him.
Data set 'a'

Note that between 0 and 30hr for all intervals of time not represented in this table there is a zero assigned as the instantaneous flow rate.
Data set 'b'

Note that between 0 and 50hr for all intervals of time not represented in this table there is a zero assigned as the instantaneous flow rate.
Solution
Each set of data is plotted, Flow Rate vs. Time. The two graphs show crackling trends.
Graph for data set 'a'
Graph for data set 'b'

Upon observing this, the general noise trend is crackling. To provide a final proof we apply the universality relation (Equation 2.5.1 ).
We first determine the critical exponent, a , for data set a.
We see that;
< S >= 1.54
S0 = 2
this value can be the first value attained or an average of the first 2 or 3 values.
T = 30
Solving for the critical exponent using the relation given above we get;
a = −0.08
Once again we carry out the same calculation for data set 'b'.
< S >= 2.12
S0 = 3
this value can be the first value attained or an average of the first 2 or 3 values
T = 50
Solving for the critical exponent using the relation given above we get;
a = −0.08
The similarity of the two critical exponent gives further proof that the data the instantaneous flow rate of the coolant over a period of
time will show crackling.
 EXAMPLE 2.5.2: COLORS OF NOISE
A chemical engineer is reading flow rates from flow meter. Every 0.1 day for 8 days a reading was taken and the data is given here:
Colors of Noise Example. The data displays the fluctuations from the set flow rate of 3000 liters per hour at a wastewater treatment
plant. The specifications for the plant say that the max flow rate is 6000 liters per hour, or the pipes will burst. Also, the flow rate
cannot fall below 200 liters per hour, or the system will automatically shut down.
The chemical engineer notices that there were some readings close to these limits. It is the chemical engineer’s job to determine if the
readings are accurate flow rates or if there is an error with the flow meter. By characterizing the type of noise, the chemical engineer
can determine the source of the noise and so take appropriate preventative measures. What type of noise is present and what protective
measures should be taken?
Solution
1) Plot data. From the data presented, a flow rate vs. time chart was graphed to gauge degree of fluctuation in the flow rate.
2) Calculate the power spectral density data. Using the simplified integral derived to calculate the power spectral density, a table was
created with the resulting PSD values at varying frequencies. These frequencies were self defined and were chosen to encompass a
broad range. More detailed calculations can be found here: Colors of Noise Example under the worksheet titled “PSD CALC”.
3) Plot the power spectral density. The power spectral density for each frequency was then plotted against the frequency, creating the
power spectral density plot shown below.
4) Characterize the noise. To determine the β value for the data, a linear trend line was taken on the data. This trend line can be seen
in the power spectral density plot above. The slope of this trend line is the β value for the data. In this case, the β=0.023. Since this
value is not that of white noise (β=0) nor that of pink noise (β=1), we can say that this noise is somewhere between white and pink
noise.
5) Determine the source of noise.
There are two possible major sources of noise. They are the liquid motion in the pipe and noise in the flow meter caused by itself or
outside sources. Since it is found earlier that β=0.023, the source of noise is probably from the flow meter. It is not from the motion of
liquid in the pipe because liquid motion tends to produce brown noise (β=2).
6) Protective measures that can be taken. Knowing this, a correction can be made on the calculations to determine the actual flow
rate. A full step by step solution can also be found here: Colors of Noise Example under the worksheet titled “Solution”.
2.5.5: SUMMARY
Noise is all around us and is an unavoidable characteristic of monitoring signals. When noise is frequency based, like the colors of noise, it
can be characterized and modeled using mathematical tools. When noise is not frequency based, like pops, snaps, and crackles, it can be
recognized, but ways to model it are still being devised. Characterizing and modeling noise is important for chemical engineers to
accurately handle their data, control processes, and predict future trends.
2.5.6: REFERENCES
Kosko, Bart. “White Noise Ain’t So White”, Noise. ISBN 0670034959
Ziemer, Rodger E. Elements of Engineering Probability and Statistics, New Jersey: Prentice Hall. ISBN 0024316202
Papoulis, Athanasios. Probability, Random Variables, and Stochastic Processes, New York: McGraw – Hill Book Company. ISBN
0071199810
Peebles, Peyton Z. Jr. Probability, Random Variables, and Random Signal Principles, New York: McGraw – Hill, Inc. ISBN 0071181814
Sethna, James P. Crackling Noise, Nature 410, 242-250 (2001). Also available here.
en.Wikipedia.org

Authors: Danesh Deonarain, Georgina Mang, Teresa Misiti, Carolyn Ehrenberger
Date Presented: September 08, 2006 Date Revised: September 19, 2006
This page titled 2.5: Noise Modeling - White, Pink, and Brown Noise, Pops and Crackles is shared under a CC BY 3.0 license and was authored, remixed,
and/or curated by Peter Woolf et al. via source content that was edited to the style and standards of the LibreTexts platform; a detailed edit history is
available upon request.
2.6: NUMERICAL ODE SOLVING IN EXCEL- EULER’S METHOD, RUNGE KUTTA,
DEAD TIME IN ODE SOLVING
Authors: (Presented: 9/8/06 /Date Revised: 9/19/06) Aaron Bennick, Bradley Anderson, Michael Salciccioli
Stewards: (9/5/07) Sean Gant, Jay Lee, Lance Dehne, Kelly Martin
2.6.1: INTRODUCTION
This article focuses on the modeling of ordinary differential equations (ODEs) of the form:
dy
= f(x, y)
dx
In creating a model, a new value y i+1is generated using the old or initial value yi, the slope estimate φ, and the step size h. This general
formula can be applied in a stepwise fashion to model the solution. All stepwise models will take the following general form:
yi+1 = yi + ϕh
The modeling methods discussed in this article are Euler’s method and the Runge-Kutta methods. The difference between the two methods
is the way in which the slope φ is estimated.
2.6.1.1: CHOOSING THE RIGHT MODEL AND STEP SIZE

The proper numerical modeling method heavily depends on the situation, the available resources, and the desired accuracy of the result. If
only a quick estimate of a differential equation is required, the Euler method may provide the simplest solution. If much higher accuracy is
required, a fifth-order Runge-Kutta method may be used. Engineers today, with the aid of computers and excel, should be capable of
quickly and accurately estimating the solution to ODEs using higher-order Runge-Kutta methods.
In all numerical models, as the step size is decreased, the accuracy of the model is increased. The tradeoff here is that smaller step sizes
require more computation and therefore increase the amount of time to obtain a solution. A balance between desired accuracy and time
required for producing an answer can be achieved by selecting an appropriate step size. One suggested algorithm for selecting a suitable
step size is to produce models using two different methods (possibly a second and third order Runge-Kutta). The steps size can then be
systematically cut in half until the difference between both models is acceptably small (effectively creating an error tolerance).
2.6.2: EULER'S METHOD

Euler's method is a simple one-step method used for solving ODEs. In Euler’s method, the slope, φ, is estimated in the most basic manner
by using the first derivative at xi. This gives a direct estimate, and Euler’s method takes the form of
yi+1 = yi + f (xi , yi ) h
dy
For demonstration, we will use the basic differential equation dx
= 3x
2
+ 2x + 1 with the initial condition y(0) = 1. If a step size, h, is
taken to be 0.5 over the interval 0 to 2, the solutions can be calculated as follows:
x y_Euler y_Actual Error
y(0) 1 1 0
y(0.5) 1 + [3(0)^2 + 2(0) + 1](0.5) = 1.5 1.875 0.375
y(1) 1.5 + [3(0.5)^2 + 2(0.5) + 1](0.5) = 2.875 4 1.125
y(1.5) 2.875 + [3(1)^2 + 2(1) + 1](0.5) = 5.875 8.125 3.25
y(2.0) 5.875 + [3(1.5)^2 + 2(1.5) + 1](0.5) = 11.25 15 3.75
The y_actual values in this table were calculated by directly integrating the differential equation, giving the exact solution as:
2.6.2.1: USE OF MICROSOFT EXCEL

Calculating an ODE solution by hand with Euler's method can be a very tedious process. Fortunately, this process is greatly simplified
through the use of Microsoft Excel. Creating a spreadsheet similar to the one above, where the x values are specified and the y_Euler values
are recursively calculated from the previous value, makes the calculation rather simple. Any step size and interval can be used.
2.6.2.2: DRAWBACKS OF EULER'S METHOD

There are several drawbacks to using Euler’s method for solving ODEs that must be kept in mind. First of all, this method does not work
well on stiff ODEs. A stiff ODE is a differential equation whose solutions are numerically unstable when solved with certain numerical
methods. For stiff equations - which are frequently encountered in modeling chemical kinetics - explicit methods like Euler's are usually
quite inefficient because the region of stability is so small that the step size must be extremely small to get any accuracy. In a case like this,
an implicit method, such as the backwards Euler method, yields a more accurate solution. These implicit methods require more work per
step, but the stability region is larger. This allows for a larger step size, making the overall process more efficient than an explicit method. A
second drawback to using Euler's Method is that error is introduced into the solution. The error associated with the simple example above is
shown in the last column. This error can be seen visually in the graph below.
It can be seen that the two values are identical at the initial condition of y(0)=1, and then the error increases as the x value increases and the
error propagates through the solution to x = 2. The error can be decreased by choosing a smaller step size, which can be done quite easily in
Excel, or by opting to solve the ODE with the more accurate Runge-Kutta method.
2.6.3: RUNGE-KUTTA METHODS

The Runge-Kutta method for modeling differential equations builds upon the Euler method to achieve a greater accuracy. Multiple
derivative estimates are made and, depending on the specific form of the model, are combined in a weighted average over the step interval.
The order of the Runge-Kutta method can range from second to higher, depending on the amount of derivative estimates made. The second-
order Runge-Kutta method labeled Heun's technique estimates derivatives by averaging endpoint measurements of the step size along a
function. This averaged value is used as the slope estimate for xi + 1. Third and higher power Runge-Kutta methods make mid-point
derivative estimations, and deliver a weighted average for the end point derivative at xi + 1. As the Runge-Kutta order increases, so does the
accuracy of the model.
The general form of the Runge-Kutta method is
yi+1 = yi + ϕ (xi , yi , h) h
Where φ(x i, yi , h) now represents a weighted average slope over the interval h .
ϕ (xi , yi , h) = a1 k1 + a2 k2 + … + an Kn
where the a's are constants and the k's are
The constants a, p, and q are solved for with the use of Taylor series expansions once n is specified (see bottom of page for derivation). The
resulting set of equations have one or more degrees of freedom. This means that for every order of Runge-Kutta method, there is a family of
methods. Below are some of the more common Runge-Kutta choices.
2.6.3.1: SECOND-ORDER RUNGE-KUTTA METHODS (N = 2)

Every second order method described here will produce exactly the same result if the modeled differential equation is constant, linear, or
quadratic. Because this is typically not the case, and the differential equation is often more complicated, one method may be more suitable
than another.
2.6.3.1.1: HEUN'S TECHNIQUE
The second-order Runge-Kutta method with one iteration of the slope estimate , also known as Heun's technique, sets the constants
1
a1 = a2 =
2
and
p1 = q11 = 1.
Huen determined that defining a and a as 1/2 will take the average of the slopes of the tangent lines at either end of the desired interval,
1 2
accounting for the concavity of the function, creating a more accurate result. When substituted into the general form, we find
1 1
yi+1 = yi + ( k1 + k2 ) h
2 2
with k1 and k2 defined as

k1 = f (xi , yi )
k2 = f (xi + h, yi + hk1 )
dy
For demonstration of this second-order Runge-Kutta method, we will use the same basic differential equation dx
2
= 3x + 2x + 1 with the
initial condition y(0) = 1. If a step size, h, is taken to be 0.5 over the interval 0 to 2, the solutions can be calculated as follows:
Figure 1: second-order Runge-Kutta
x k_1 k_2 y_Heun y_Actual Error
y(0) 3(0)^2 + 2(0) + 1 = 1 3(0.5)^2 + 2(0.5) + 1 = 2.75 1 1 0
y(0.5) 3(0.5)^2 + 2(0.5) + 1 = 2.75 3(1)^2 + 2(1) + 1 = 6 1 + [0.5(1) + 0.5(2.75)](0.5) = 1.9375 1.875 -0.0625
y(1) 3(1)^2 + 2(1) + 1 = 6 3(1.5)^2 + 2(1.5) + 1 = 10.75 2.375 + [0.5(2.75) + 0.5(6)](0.5) = 4.125 4 -0.125
y(1.5) 3(1.5)^2 + 2(1.5) + 1 = 10.75 3(2)^2 + 2(2) + 1 = 17 5.375 + [0.5(6) + 0.5(10.75)](0.5) = 8.3125 8.125 -0.1875
y(2.0) 3(2)^2 + 2(2) + 1 = 17 3(2.5)^2 + 2(2.5) + 1 = 24.75 10.75 + [0.5(10.75) + 0.5(17)](0.5) = 15.25 15 -0.25
When compared to the Euler method demonstration above, it can be seen that the second-order Runge-Kutta Heun's Technique requires a
significant increase in effort in order to produce values, but also produces a significant reduction in error. Following Runge-Kutta methods
can be worked through a similar manner, adding columns for additional k values. Below is a graphical description of how slope is estimated
using both Euler's method, and Heun's technique. Observe the increase in accuracy when an average slope across an interval of 0.5 is used
instead of just an initial estimate.
2.6.3.1.2: IMPROVED POLYGON METHOD
Using the improved polygon method, a2 is taken to be 1, a1 as 0, and therefore . The general form then becomes
Ralston's Method
The Ralston method takes a2 to be . Therefore and . It has been determined by Ralston (1962) and Ralston and
Rabinowitz (1978) that defining a2 as will minimize truncation error in second-order Runge-Kutta methods. The general form becomes
2.6.3.2: THIRD-ORDER RUNGE-KUTTA METHODS (N = 3)

The third-order Runge-Kutta methods, when derived, produce a family of equations to solve for constants with two degrees of freedom.
This means an even more variable family of third-order Runge-Kutta methods can be produced. A commonly used general third-order form
is
with
2.6.3.3: FOURTH-ORDER RUNGE-KUTTA METHODS (N = 4)

The family of fourth-order Runge-Kutta methods have three degrees of freedmon and therefore infinite variability just as the second and
third order methods do. The fourth-order versions are most favored among all the Runge-Kutta methods. What is know as the classical
fourth-order Runge-Kutta method is
with
2.6.3.4: FIFTH-ORDER RUNGE-KUTTA METHODS (N = 5)

Where very accurate results are required, the fifth-order Runge-Kutta Butcher's (1964) fifth-order RK method should be employed:
with
The integration of k's within other k values suggests the use of a spreadsheet. As with all Runge-Kutta methods, the calculation of values for
the fifth-order version would be greatly assisted through the use of Microsoft Excel.
2.6.4: INTERACTIVE MODEL COMPARISON IN EXCEL

To give a better understanding of the impact between different Runge-Kutta methods, as well as the impact of step size, the interactive excel
sheet below will allow you to enter the step size h into a set of models and observe how the models contour to match an example differential
equation solution.
Observe the relationship between model type, step size, and relative error.
ODE model comparison interactive spreadsheet
2.6.5: DEAD TIME

Dead time, or delay differential equation, occurs when there is a delay or lag in the process of a real life function that is being modeled. This
means that the numerical model is not accurate until the delay is over. This concept can come into play for the start up of a reaction process.
For example, at the start of a reaction in a CSTR (Continuous Stirred Tank Reactor), there will be reagents at the top of the reactor that have
started the reaction, but it will take a given time for these reactant/products to be discharged from the reactor. If one was modeling the
concentration of reagents vs. time, time t=0 would have started when the tank was filled, but the concentrations being read would not follow
a standard model equation until the residence time was completed and the reactor was in continuous operational mode. The dead time is the
time it would take for the readings to start meeting a theoretical equation, or the time it takes for the reactor to be cleared once of the
original reagents. Dead time can be determined experimentally and then inserted into modeling equations.
The solution of such a delay differential equation becomes problematic because in order to solve the equation, information from past times
(during the delay) is needed in addition to the current time. For example if to is the lag time for the given scenario, then the value to use in
the ODE becomes (t-to) instead of t. When modeling this in excel, (x-to) is substituted in for the x value. With this substitution, Euler’s
method can be used again in the same way to approximate the solution. This substitution will be carried into the differential equation so that
the new dead time solution can be approximated using Euler's method. To model the reactor before the initial deadtime is completed, piece-
wise functions are often used. This way for the dead time, a given model is used not characteristic of the reactor at normal operating
conditions, then once the dead time is completed the modeling equation is taken into effect.
The link below will help to show how to include dead time in a numerical method approximation such as Euler's method. As seen in the
excel file, the dead time that is specified by the user in the yellow box will change the delay in the model. The more dead time, the further
shifted from the theoretical equation the new model is. To take dead time into account in excel, the x value is simply substituted out for (x-t)
where t is equal to the dead time. If you look at the equations entered in the Y cells, you will see that the x value inserted into the differential
equation is (x-t), where t is the user specified dead time. It can be seen through this example spreadsheet that the effect of dead time is a
simple horizontal shift in the model equation.
Dead Time Interactive Spreadsheet
2.6.6: ERROR
There are two types of error associated with solving ODEs using stepwise approximation methods in Excel. These errors are also present
using other methods or computer programs to solve ODEs. The first, discretization, is the result of the estimated y value that is inherent of
using a numerical method to approximate a solution. Discretization errors, also called truncation, occur proportionately over a single step
size. Truncation error will propagate over extended results because the approximation from previous steps is used to approximate the next.
In essence, a copy of a copy is being made. The accuracy of the next point is a direct result of the accuracy of the previous. Just as the
quality of a second copy is dependant on the quality of the first. This error can be reduced by reducing the step size.
Please see ODE model comparison interactive spreadsheet to better learn how step sizes can influence error.
The second types of errors are rounding errors. These errors are dependent on the computer’s capacity for retaining significant digits. The
more significant digits that a computer can hold, the smaller the rounding error will be.
2.6.6.0.1: ESTIMATING ERROR IN EULER'S METHOD
To mathematically represent the error associated with Euler's method, it is first helpful to make a comparison to an infinite Taylor series
expansion of the term yi + 1. The Taylor series expansion of this term is
2 3 n
′
h ′′
h n
h
yi+1 = yi + f (xi , yi ) ℏ + f (xi , yi ) +f (xi , yi ) +…+f (xi , yi )
2 3 !n
When this expansion is compared to the general form of Euler's method it can be seen that Euler's method lacks every term beyond
. These missing terms, the difference between the Euler approximation and an infinite Taylor series (taken to be the true
solution), is the error in the Euler approximation. Mathematically respresenting the error in higher order Runge-kutta methods is done in a
similar fashion.
2.6.6.0.1: ESTIMATING AND MINIMIZING ERROR IN RUNGE KUTTA METHOD
Uniform time steps are good for some cases but not always. Sometimes we deal with problems where varying time steps makes sense.
When should you change the step size? If we have an nth order scheme and and (n+1)th order scheme, we can take the difference between
these two to be the error in the scheme, and make the step size smaller if we prefer a smaller error, or larger if we can tolerate a larger error.
This is fairly simple with Runge Kutta, because we can take a fifth order method and a fourth order method using the same k's. Only a little
extra work at each step.
Another way of estimating error in the Runge-Kutta method is to reverse directions at each step of the advancing solution and recompute the
previous ordinate. By considering the difference between the newly computed previous ordinate and the originally computed value, you can
determine an estimate for the truncation error incurred in advancing the solution over that step. This is a bit more tedious, but does give a
good estimate of truncation error.
2.6.7: EULER'S METHOD FOR SYSTEMS OF ODES

In chemical engineering and other related fields, having a method for solving a differential equation is simply not enough. Many real world
problems require simultaneously solving systems of ODEs. For problems like these, any of the numerical methods described in this article
will still work. The only difference is that for n ODEs, n initial values of y are needed for the initial x value. Then the numerical method of
choice is just applied to every equation in each step before going on to the next. This further complicates the step-by-step problem solving
methodology, and would require the use of Excel in nearly every application. The ODE solved with Euler's method as an example before is
now expanded to include a system of two ODEs below:
dy1
2
= 3x + 2x + 1, y1 (0) = 1
dx
dy2
= 4y1 + x, y2 (0) = 2
dx
Step size is again 0.5, over an interval 0-2.

x y_1Euler y_2Euler y_1Actual y_2Actual y_1 Error y_2 Error
y(0) 1 2 1 2 0 0
y(0.5) 1 + [3(0.5)^2 + 2(0.5) + 1](0.5) = 2.375 2 + [4(1) + 0.5](0.5) = 4.25 1.875 5.875 -0.5 1.625
y(1) 2.375 + [3(1)^2 + 2(1) + 1](0.5) = 5.375 4.25 + [4(2.375) + 1](0.5) = 9.5 4 18.5 -1.375 9
y(1.5) 5.375 + [3(1.5)^2 + 2(1.5) + 1](0.5) = 10.75 9.5 + [4(5.375) + 1.5](0.5) = 21 8.125 51.875 -2.625 30.875
y(2.0) 10.75 + [3(2)^2 + 2(2) + 1](0.5) = 19.25 21 + [4(10.75) + 2](0.5) = 43.5 15 124 -4.25 80.5
It should be observed that there is more error in the Euler approximation of the second ODE solution. This is because the equation also has
y1 in it. So there is the error introduced by using the Euler approximation to solve the 2nd ODE, as well as the error from the Euler
approximation used to find y1 in the 1st ODE in the same step!
Exact solutions were again obtained from directly integrating the ODEs: and
2.6.7.1: EXAMPLE 1: REACTOR DESIGN

The elementary liquid-phase reaction A --> B is to be carried out in an isothermal, isobaric PFR at 30 degrees C. The feed enters at a
concentration of 0.25 mol/L and at a rate of 3 mol/min. The reaction constant is known experimentally to be 0.01 min-1 at this temperature.
You have been given the task of building a reactor that will be used to carry out this reaction. Using Euler’s method with a step size of 0.05,
determine how large the reactor must be if a conversion of 80% is desired.
Simplified Design Equation:
dV FA0
=
dX k CA0 (1 − X)
Lumping the given flow, concentration, and reaction constant together gives:
dV 1
= 1200
dX 1− X
Since no volume is required for a conversion of zero, the initial condition needed is V(0)=0.
Now the information simply has to be entered into Excel. A column for the conversion, X, going from 0 to 0.8 in 0.05 increments is used for
the step size, and the first value, V(0), is known to be zero. To get the volume, simply add the previous volume to the constants multiplied
by the step size and 1/(1-X), or:
1
Vi+1 = Vi + 1200 ∗ 0.05 ∗
1− X
Copying this formula down the column to the final conversion value of 0.8, gives the results shown in the table below:
X V
V(0) 0
V(0.05) 63
V(0.1) 130
V(0.15) 200
V(0.2) 275
V(0.25) 355
V(0.3) 441
V(0.35) 533
V(0.4) 633
V(0.45) 743
V(0.5) 863
V(0.55) 996
V(0.6) 1146
V(0.65) 1317
V(0.7) 1517
V(0.75) 1757
V(0.8) 2057
The final reactor volume obtained is approximately 2057 L. This compares reasonably well to the exact value of 1931 L. The Excel file
used to obtain this solution, along with the exact solution, can be downloaded below.
Example 1
2.6.8: SECOND-ORDER RUNGE-KUTTA DERIVATION

The following example will take you step by step through the derivation of the second-order Runge-Kutta methods. Setting n = 2 results in a
general form of
The constants in the general form must be defined. To do this we will employ a second-order Taylor series expansion for yi + 1 in terms of yi
and . This Taylor series is
2
′
h
yi+1 = yi + f (xi , yi ) h + f (xi , yi )
2
Expanding with the chain rule, and substituting it back into the previous Taylor series expansion gives
2
∂f ∂f ∂y h
yi+1 = yi + f (xi , yi ) h + ( + )
∂x ∂y ∂x 2
The next step is to apply a Taylor series expansion to the k2 equation. The applied Taylor series expansion rule is
∂g ∂g
g(x + r, y + s) = g(x, y) + r +s +…
∂x ∂y
and the resulting equation is

∂f ∂f
2
f (xi + p1 h, yi + q11 k1 h) = f (xi , yi ) + p1 h + q11 k1 h + O (h )
∂x ∂y
where O(h2) is a measure of the truncation error between model and true solution.
When this Taylor series expansion result of the k2 equation, along with the k1 equation, is substituted into the general, and a grouping of like
terms is performed, the following results:
∣ ∂f
∣a q ∂f
∣ 2 3
yi+1 = yi [a1 f (xi , yi ) |a2 f (xi , yi )] h [a3 p ∣ s 11 f (xi , yi ) ∂ y ∣ h |O ( h )
∣ 1
Setting this result equal to the substituted general form will allow us to determine that, for these two equations to be equivalent, the
following relationships between constants must be true.
It should be noticed that these three equations, relating necessary constants, have four unknowns. This means that to solve for three
constants, one must first be chosen. This results in a family of possible second-order Runge-Kutta methods.

Example of usage of Euler's Method:
video.google.com/googleplayer.swf?docId=1095449792523736442
Narrated example of using the Runge-Kutta Method:
video.google.com/googleplayer.swf?docId=-2281777106160743750
Unnarrated example of Using the Runge-Kutta Method:
File:Numerical Solving in Excel, Unnarrated.ppt
2.6.10: REFERENCES
"Delay Differential Equation", Wolfram MathWorld, Online: September 8, 2006. Available
http://mathworld.wolfram.com/DelayDifferentialEquation.html
Chapra, Steven C. and Canale, Raymond P. "Numerical Methods for Engineers", New York: McGraw-Hill.
Franklin, Gene F. et al. "Feedback Control of Dynamic Systems", Addison-Wesley Publishing Company.
R. England. "Error Estimates for Runge-Kutta Type Solutions to Systems of Ordinary Differential Equations", Research and
Development Department, Pressed Steel Fisher Ltd., Cowley, Oxford, UK. October 1968.
Call, Dickson H. and Reeves, Roy F. "Error Estimation in Runge Kutta Procedures", ACM, New York, NY, USA. September 1958.
This page titled 2.6: Numerical ODE solving in Excel- Euler’s method, Runge Kutta, Dead time in ODE solving is shared under a CC BY 3.0 license and
was authored, remixed, and/or curated by Peter Woolf et al. via source content that was edited to the style and standards of the LibreTexts platform; a
detailed edit history is available upon request.
2.7: SOLVING ODES WITH MATHEMATICA- HOW TO FIND NUMERICAL AND
ANALYTICAL SOLUTIONS TO ODES WITH MATHEMATICA
Authors: Matthew Baumgartner, Olyvia Dean, Viral Patel, Joel Schweitzer, and Eric Van Beek
Stewards: Brian Hickner, Lennard Gan, Addison Heather, Monique Hutcherson
Date Released: September 6, 2006 /Date Revised: September 8, 2007
2.7.1: INTRODUCTION
Mathematica is an advanced mathematics solution program created by Wolfram Research, Inc. One of the most powerful software packages
of its kind, Mathematica is able to easily and conveniently solve complicated mathematical formulae, including differential equations. This
article focuses on the modeling of first and higher order Ordinary Differential Equations (ODE) in the following forms:
dy
= f(x, y) (Basic ODE)
dx
z
dy
= f(x, y) (Higher Order ODE)
dz x
Like all software, Mathematica uses a specific language in order to perform calculations. The names of all functions must be capitalized-
this applies to both mathematical functions (such as Sin and Cos) and built-in functions (such as Plot and DSolve). For this reason it is
common for users to write their own functions in minimized letters. This decreases the chance of overriding or redefining a Mathematica
function. Square brackets always follow function calls, and the function's parameters are always enclosed in curved brackets. For example,
if a user wished to plot sin(x) from x = 0 to x = 1, they would type: Plot[Sin[x],{x,0,1}]. The user must type "Shift"+"Enter" to input the
function. Typing only "Enter" will add a line to the formula. . For full PC keyboards, the “Enter” key on the far right is the equivalent of
“Shift”+”Enter”. If the user does not want to see the output line, it can be suppressed by typing a semi-colon at the end of the expression.
Mathematica also recognizes common mathematical constants such as pi (Pi), Euler's constant (E), and imaginary numbers (I). Note that
these constants must also be capitalized when entered.
Mathematica features two functions for solving ODEs: DSolve and NDSolve. DSolve is used when the user wishes to find the general
function or functions which solve the differential equation, and NDSolve is used when the user has an initial condition. The prompts for the
functions are quite similar. Note: Two equal signs must be used to denote equality in an equation. Using one equal sign assigns a value to a
variable.
Example:
f(x) = 5x2 + 7 This input creates a variable named f(x) that has a value of 5x2 + 7.
f(x) = = 5x2 + 7 This input creates a function f(x) which is defined by 5x2 + 7.
Mathematica will sometimes store user created functions or variables into its cache. This causes some (potentially) correct commands to
output errors. It is recommended to quit the local kernel and reinitialize the commands after any changes to your code.
Also, check that every command in the code is being executed. Clearing all the outputs in the cell may be helpful to figure out which
commands had not been executed. This can be done by going to the “Cell” option found on the top and choosing the “Delete All Output”
Example:
OUTPUT IS NOT CLEARED
CLEARING THE OUTPUTS
CODE WITH CLEARED OUTPUTS
As seen in the examples, it is easier to troubleshoot and debug your program when it looks less confusing. Clearing the extra outputs helps
you focus on just the code that you have written.
2.7.2: FIRST ORDER ODES

Notation: In the following examples, eqn represents the ODE, y represents the function being solved for, i represents the initial condition,
and x and t are independent variables.
2.7.2.1: ODES WITH INITIAL CONDITIONS

NDSolve will not display a numerical value, but rather an “Interpolating Function” which can be displayed graphically. The easiest way to
display this graph is to assign the solution to a variable (called Solution in the example) in the input line, and then use the “Plot” function to
display it.
Input for One ODE: Solution = NDSolve[{eqn,y[0] == i},y,{x,xmin,xmax}]
Input for Multiple ODEs: Solution = NDSolve[{eqn1,eqn2,...,y1[0] == i1,y2[0] == i2,...},{y1,y2,...},{x,xmin,xmax}]
Input for Partial Differential Equations: Solution = NDSolve[{eqn,y[0] == i},y,{x,xmin,xmax},{t,tmin,tmax}]
After Mathematica solves the ODE, plot the solution by typing: Plot[Evaluate[y[x] /. Solution],{x,xmin,xmax}], where “/.” is a notation
used by Mathematica meaning “following the rule of”. It basically recalls the function stored as the variable in the previous input line.
Depending on the needs of the user, the functions “ParametricPlot” or “Plot3D” may also be useful. The notation for these functions is the
same as “Plot”.
Here is a simple example of what you would type into Mathematica:
Solution = NDSolve[{y'[x]==y[x]*Cos[x+y[x]],y[0]==1},y,{x,0,30}]
Mathematica will output:
Output[1]= {{y->InterpolatingFunction[ { { 0.,30.} },<>] } }
To plot this function you would type:
Plot[Evaluate[y[x]/.Solution],{x,0,30}]
Note: Remember to type "Shift"+"Enter" to input the function
2.7.2.2: ODES WITHOUT INITIAL CONDITIONS

Input for One ODE: DSolve[eqn,y,x]
Input for Multiple ODEs: DSolve[{eqn1,eqn2,…},{y1,y2,…},x]
Input for a Partial Differential Equation: DSolve[eqn,y,{x1,x2,…}]
Solution = DSolve[y'[x] ==4*x-2*x*y[x],y[x],x]
Output[1]=
ParseError: invalid Primary (click for details)
2.7.3: SECOND ORDER AND HIGHER ODES

Input for Higher Order ODEs: Solution = NDSolve[{eqn,y[0] == i1, y' [0] == i2,…},y,{ x,xmin,xmax}]
Plotting the solution: Plot[Evaluate[{y[x], y' [x],...} /. Solution],{x,xmin,xmax}]
Solution = NDSolve[{y[x] +Sin[y[x]]+y[x]==0,y[0]==1,y'[0]==0,y,{x,0,30}]
Output[1]= {{y->InterpolatingFunction[ { {0.,30.} },<>] } }
To plot this function you type:
Plot[Evaluate[{y[x],y'[x],y[x]}/.Solution],{x,0,30}]
Mathematica should output a graph.
2.7.4: ALGORITHMS USED BY MATHEMATICA

Mathematica uses two main algorithms in order to determine the solution to a differential equation. These algorithms are the Adams method
and the Gear method.
The Adams and Gear methods are forms of linear multistep methods. An example of these would be the following:
In the example above, h denotes the step size and the coefficients are determined by the method used. Multistep methods are
expansions of more familiar single-step methods used to solve differentials (i.e. Euler, Runge-Kutta, Taylor). Each of these methods
requires an initial point in order to calculate the following point. Similarily, multistep methods also require initial points in order to solve the
ODE. The number of initial points required depends on which method is used to solve the ODE. Multistep methods typically produce less
error than the single-step methods because of multiple initial points.
In order to determine what method to use one must first find the stiffness of the function.
2.7.4.1: ADAMS METHOD

Euler, Taylor and Runge-Kutta methods used points close to the solution value to evaluate derivative functions. The Adams-Bashforth
method looks at the derivative at old solution values and uses interpolation ideas along with the current solution and derivative to estimate
the new solution [4]. In order to solve an ODE using this method, f(t,y) must be continuous and satisfy Lipschitz condition for the y-variable
which states [5]:
for all | h | < ε where B and β are independent of h, β > 0 and α is an upper bound for all β for which a finite B exists
This is a basic form of the Adams-Bashforth method. Note that two initial points are required for this method.
There is another Adams method that requires three initial points. The method is solved the same way, however the equation varies a little bit
and is referred to as the Adams-Moulton Method.
The coefficients/constraints, β can be solved for using knowledge of ODE's and other math tools. It was stated earlier that .
We can let f(t,y) = λy, therefore . We can also let if there is a constant step size and σ represents a polynomial.
Through substitution we find [6]:
We can expand the quadratic using another math identity and ultimately solve for constraints β1 and β2. Another method for solving for
coefficients β1,β2 is mentioned below:
In order to find the coefficient βj one must first use polynomial interpolation to find the polynomial p of degree s − 1 such that:
for
From this the Lagrange formula for polynomial interpolation yields
Now the polynomial p is a locally good approximation of the right-hand side of the differential equation y' = f(t,y) that is to be solved. Now
we must consider the equation y' = p(t) instead. This equation can be solved exactly by simply taking the integral of p.
The Adams–Bashforth method arises when the formula for p is substituted. The coefficients bj turn out to be
s−j 1 s
(−1)
bet aj = ∫ ∏(u + i) du.
(j − 1)!(s − j)! 0 i=1
i≠j
The Adams-Bashforth method is typically used for Linear and Non-liner ODE's with dense systems.
2.7.4.2: GEAR METHOD

Taking
bet ai = 0
for i ≥ 1
The Gear method, also known as the backward differentiation formulae (BDF, a.k.a. Gear’s formulae) is another multi-step method but is
generally used for multiple equations. In order to use the gear method your function must have a stiffness greater than 500, meaning that the
function is stiff.
As a particular cases, taking for i ≥ 2 and optimizing the remaining coefficients to maximize the accuracy of the resulting scheme
recovers the Implicit Euler method. Taking for i ≥ 3 gives
We now focus on this case in particular. Applying this method to the scalar model problem and assuming constant h and a solution of the
form , we find the following quadratic equation for
the two roots of which are given by
− −−−−−− − −−−−−−
− α1 + √ γ(1 + ϵ) − α1 − √ γ(1 + ϵ)
sigma = or
2(1 − Σ) 2(1 − Σ)
2
gamma = α1 − 4α2
Applying the identities, we may expand both roots in terms of powers of h. By our assumed form of the solution, it follows that
.The leading-order term in the expansion in h of (a “spurious root”) is proportional to h. For small h, quickly
decays to zero, and thusmay be neglected. The leading-order terms in the expansion in h of (the “physical root”) resemble the Taylor-
series expansion of the exact solution over a single timestep.
Matching coefficients with the expansion of the exact solution as indicated by underbraces in the above
expression, and applying the definition , we arrive at three equations for to achieve the highest order of
accuracy possible with this form. It is easily verified that satisfy these three equations. The leading-
order error term of this method is proportional to . Thus, over a single timestep, the scheme is “locally third-order accurate”; more
significantly, over a fixed time interval [0,T], the scheme is globally second-order accurate. The resulting method,
is thus referred to as BDF2, and may be viewed as a implicit alternative to AM3 that, for the same number of steps into the past, p =
max(m,n), has reduced order of accuracy but greatly improved domain of stability. Higher-order BDFs may be derived in an analogous
fashion; BDF3, BDF4, BDF5, and BDF6 in particular are found to have excellent stability properties as compared with their AM
counterparts with the same number of steps.

You have been assigned a non-isothermal reactor with the following material and enthalpy balances.
2 2
fracdca dt = ca (t ) − cb (t) − T (t )
Species a Material balance

3 2
fracdcb dt = cb (t ) − ca (t) − T (t )
Species b Material balance

2
fracdT dt = ca (t) − cb (t) − T (t )
Enthalpy balance
When is the concentration of species a and species b equal? The initial conditions for the reactor are as follows:
T0 = 1
ca0 = 2.0
cb0 = 1.8
tf = 0.2
s = NDSolve[{x'[t] == x[t]^2 - y[t] - z[t]^2, y'[t] == y[t]^3 - x[t] - z[t]^2, z'[t] == x[t] - y[t] - z[t]^2, x[0] == 2, y[0] == 1.8, z[0] == 1},
{x, y, z}, {t, 0.2}]
Plot[Evaluate[{x[t],y[t]} /. s], {t, 0, 0.2}]
You have been asked to study a semi-batch reactor for the reaction . The rate law is , where k = 2.7. Other
parameters are: initial volume = 5, volumetric flow rate = 0.05, initial concentration of A = 2500. Use Mathematica to create a Conversion
versus Time plot for 100 time units.
1. Mole Balance:
ra V
fracdXdt = −
NA0
2. Rate Law:
3. Stoichiometry:
4. Combine:
If the semi-batch reactor runs for 100 time units, the conversion is about 0.8. This can be seen in the plot below or by running the two
commands that follow into Mathematica.
s = NDSolve[{y'[t] == 2.7 (1 - y[t]) (5 + .05 t)/(2500 y[t]), y[0] == 0.0001}, y, {t, 0, 100}];
Plot[Evaluate[y[t] /. s], {t, 0, 100}, PlotRange -> All]
This is the Mathematica notebook file for the example: Media:MathematicaEx2.nb

Note: The parameters used in this problem are fabricated and are intended to illustrate the use of Mathematica in solving ODEs.
2.7.7: FORMATTING PLOTS IN MATHEMATICA

Here are some useful tips to formatting plots in Mathematica:
Label Plot: e.g. Plot[Sin[x], {x,0,2Pi}, PlotLabel->"Volume of Tank vs. Time"]
Label Axes: e.g. Plot[Sin[x], {x,0,2Pi}, AxesLabel->{"time(min)","volume(L)"}]
Color Plot: e.g. Plot[Sin[x], {x,0,2Pi}, PlotStyle->Red]. You can also make the plot thick, dashed, etc. e.g. Plot[Sin[x], {x,0,2Pi},
PlotStyle->{Red,Thick,Dashed}]
Insert Legend: e.g. Needs["PlotLegends`"] <--MUST insert this BEFORE your plot command. Plot[{Sin[x],Cos[x]}, {x,0,2Pi},
PlotLegend->{"sine", "cosine"}]

Solving ODEs with Mathematica Video
slides for this talk
2.7.9: ADDITIONAL TIPS AND TRICKS FOR TROUBLESHOOTING IN MATHEMATICA

Mathematica is a powerful computing tool, however the syntax can be a bit difficult to understand. Here are some notes for troubleshooting
in Mathematica.
1. Check to make sure that your variable names and signs are consistent.
Ex) Make sure you use xI everywhere instead of xI and x1 or xl.
Ex) Functions, including the ones you create, are usually followed by brackets such as Sin[x] or y[x]. However, brackets are not
necessary when you are solving for a function in a set of differential equations such as NDSolve[eqns, {y}, {x, 0, 50}];
Ex) Check to see if your parentheses are aligned such that you are actually entering the function you think you're entering. Recall order
of operations and the distributive property. x*(a+b) is NOT equal (x*a) + b. This seems simple, but often gets overlooked when dealing
with several lines of code.
2. You may find it easier to define all of your variables together at the beginning of your code. This makes it easy to go back and change the
values assigned to each variable if trying to determine the impact one variable has on a system of equations. For instance, say you are trying
to determine the effects the flow rates into two different tanks (F1, F2) will have on the tank volumes over ten time steps. The differential
equations governing the situation are: and , where F1 = 2, F2 = 12, V1(0) = 0, V2(0) = 0.
If you write the ODEs in Mathematica by directly substituting in F1 = 2 and F2 = 12, you will have to change the ODEs each time you
change the values of F1 and F2. Below is an example of what the Mathematica equations would look like.
s = NDSolve[{V1’[t] == 8 – V1[t], V2’[t] == 12 – (1/3)*V2[t], V1[0] == 0, V2[0] == 0},{V1,V2},{t,0,10}]
Another option is to define F1 and F2 before your equations, and then call them when solving for your ODEs. This allows you to easily
substitute in new values for your variables without changing the differential equations, and reduces the odds of making simple
computational errors in your ODE in the process. The Mathematica code would look as shown below.
variables = {F1 -> 2, F2 -> 12};
s = NDSolve[{V1’[t] == 4*F1 – V1[t] , V2’[t] == F2 – (1/3)*V2[t], V1[0] == 0, V2[0] == 0} /. variables,{V1,V2},{t,0,10}]
3. Be Aware of the Kernel
The Mathematica Kernel stores all of the information about variable definitions. When you define a variable , the definition is stored there.
The kernel is automatically started when a new Mathematica session is started. You may also start it manually. To start the kernel manually
go to Evaluation -> Start Kernal -> Local. Once the kernel is started and you wish to go back and change a variable definition you must
"Quit the Kernal" before you see the change occur. In order to "Quit the Kernal" you must go to Evaluation -> Quit Kernal -> Local.
Terminating the Mathematica Kernal, erases all of the definitions previously entered. For this reason, after you "Quit the Kernal" and enter
in new definitions for your variables you must re-enter all of your code. The images below show how to "Start Kernal" and how to "Quit the
Kernel."
Example:
An example of when you would want to "Quit the Kernal" when using Mathematica for Controls is when you are finding a steady state
point using a PID controller. We start by defining the variables for vset, taui, taud, and Kc. Mathematica solves this equation using the
defined variables and shows that (x->3.9, V->10) is the steady state point. This is shown in the figure 1 below.
Figure 1 when you defined Kc even though you have not defined it here.
Figure 3 shows the output once you have "Quit the Kernel" and re-entered the Mathematica code.
Now you see that the definition for Kc has been deleted because the steady state points are in terms of Kc. To find the impact of Kc you can
use this solution.
4. Define functions or formulas that are used often
If you will use a function often and it is not already defined in Mathematica by default, define your own. This can be especially helpful if
you intend repeat the same function multiple times while only changing variables. A scenario where defining your own function could save
you time and error is when comparing P-Fisher's values while keeping the margins constant. Only the variables will change but the function
will remain the same. The new function will only last as long as each session in Mathematica.
To define a function, begin by naming the function followed by an open bracket. List all the variables with an underscore after each name
and separate each variable with a comma. Close the bracket when all the variables have been listed. To let Mathematica know that this is a
newly defined function, after closing the bracket place a semi-colon and equal sign. Now define what the formula is using the variable
names. When you have finished, click shift-enter simultaneously. To test the new formula, start with naming the function and open bracket,
then list all the variables numerical values separated by comma's in the same order that was listed when defining the function. Close the
bracket after all the variables have been listed and click shift-enter simultaneously. The output will be the next line in Mathematica with the
answer.
Example:
It is also possible to use this tool to put equations in terms of unknown variables. To do this begin by defining the function the same way as
before. When plugging in the numerical values of the variables, leave the unknown variables as their variable name. The Mathematica
output will provide the answer in terms of the unknown variables.
Example:
5. "Make it pretty"
Oftentimes when people program in any language, they tend to crowd lines, almost as if they were trying to save paper. Well, since this does
not apply when you are typing on a screen, don't worry about it.
Ex) You could write:
vars = {V -> 1000, Cao -> 0.2, UA -> 20000, F -> 2000, DE1 -> 30000, k1 -> 1.0 10^13, DH1 -> 50000, rcp -> 60, R -> 1.987, To -> 500,
Tf -> 500, Tjin -> 500, Fj -> 1000, Vj -> 100, Caf -> 0.2, DH2 -> -192000, DE2 -> 30000, k2 -> 1.0 10^13};
eqns = {Ca'[t] == (Caf - Ca[t]))/V - Ca[t]*k1*Exp[-DE1/(R*T[t])], Cb'[t] == (0 - Cb[t])/V +k1*Ca[t]*Exp[-DE1/(R*T[t])] -
k2*Cb[t]*Exp[-DE2/(R*T[t])], T'[t] == (Tf - T[t])/V + (-DH1/(rcp))*k1*Ca[t]*Exp[-DE1/(R*T[t])] + (-DH2/rcp )*k2*Cb[t]*Exp[-
DE2/(R*T[t])] - (UA (T[t] - Tj[t]))/(V *rcp), Tj[t] == (Fj (Tjin - Tj[t]))/Vj + (UA (T[t] - Tj[t]))/(Vj *rcp), Ca[0] == 0.2, Cb[0] == 0, T[0] ==
500, Tj[0] == 500};
sol = NDSolve[eqns /.vars, {Ca, Cb, T, Tj}, {t, 0, 50}]; Plot[{Ca[t]} /. sol, {t, 0, 50}, AxesLabel -> {t, Ca}, PlotRange -> Full];Plot[{Cb[t]}
/. sol, {t, 0, 50}, AxesLabel -> {t, Cb},PlotRange -> Full];Plot[{T[t]} /. sol, {t, 0, 50}, AxesLabel -> {t, T}, PlotRange -> Full];Plot[{Tj[t]}
/. sol, {t, 0, 50}, AxesLabel -> {t, Tj}, PlotRange -> Full]
but it looks much better if you do this:
These thoughts will hopefully help you later on and make it easier to read your code.
6.) "Check the colors!"
Here is a list of font color that Mathmatica will output for a specific reason.
Local Variables in a certain Font Color
Local variables of Module and With in Green
1. Example:
Function arguments and pattern names in Green (Italics)
1. Example:
Variables made special by use in arguments in Turquoise
Example:
Errors and Warnings in a certain Font Color
Syntax Errors in Purple
Example:
Emphasized Syntax Errors in Dark Red with Yellow Background
Example:
Missing arguments in Bright Red
1. Example:
Excess arguments in Bright Red
1. Example:
Possible unwanted assignments in Bright Red
Example:
Unrecognized option names in Bright Red
Example:
Local Scope conflicts in Dark Red with Yellow Background
Example:
Variables that will go out of scope before being used in Bright Red
Example:
Shadowing in multiple contexts in Bright Red
Example:
Other in a certain Font Color
Comments in Light Grey
Example:
Strings in Dark Grey
Example:
Global symbols that have no value assigned in Bright Blue
Example:
2.7.10: ACCESSING MATHEMATICA FROM YOUR PERSONAL COMPUTER

If you prefer to use Mathematica from your home computer rather than one on campus, you can remotely log into a CAEN computer
through U-M’s Virtual Sites service. Follow the steps outlined below:
2. Visit http://virtualsites.umich.edu/ . The U-M WEBLOGIN page will appear. Log in using your unique name and Kerberos password.
3. The Virtual Sites page should appear on your screen. Click on the large, orange Connect Now button.
4. Under the heading What software do you need?, select Engineering (CAEN) Legacy XP. Then click on the Request Connection
button toward the bottom of the page.
5. Follow the instructions under the heading Launch Virtual Sites by downloading and opening the file in the Internet Explorer
information bar. A window will pop up on your screen titled Remote Desktop Connection. Click the Connect button.
6. You will then see the screen that appears when logging into a CAEN computer. Your unique name should already be included. Type your
Kerberos password and click OK to log in to the CAEN computer.
7. Click OK on the Notice window that pops up and your personal settings will be applied. Also click OK if a System Requirements
Wizard window appears.
8. If you scroll down using the scroll bar on the right-hand side of the screen, you will see the Start button. Click Start -> All Programs -
> Math and Numerical Methods -> Wolfram Mathematica -> Mathematica 6. Once Mathematica opens on the computer screen,
you may use it as though you were actually sitting in front of a CAEN computer on campus.
9. To log out, click Start -> Log Off. Under the Log Off Windows window, click Log Off.
2.7.11: REFERENCES
"Adams Method", Wolfram MathWorld, Online: August 5, 2007. Available http://mathworld.wolfram.com/AdamsMethod.html [1]
"Gear Predictor-Corrector Method", Davidson University, Online August 5, 2007. Available
webphysics.davidson.edu/Projects/SuFischer/node47.html [2]
Fogler, H. Scott (2006), Elements of Chemical Reaction Engineering, New Jersey: Pretince Hall PTR. ISBN 0-13-047394-4 [3]
"Adams-Bashforth", Pitt Math, Online: September 8, 2007. Available:
www.math.pitt.edu/~sussmanm/2071Spring07/lab02/index.html#AdamsBashforthMethods [4]
"Lipschitz Condition", Wolfram MathWorld, Online: September 7, 2007. Available
http://mathworld.wolfram.com/LipschitzCondition.html [5]
"Adams-Bashforth (AB) Methods", University of California-San Diego, Online: September 7, 2007. Available:
renaissance.ucsd.edu/chapters/chap11.pdf [6]
"Gears Method", University of California-San Diego, Online: September 7, 2007. Available: renaissance.ucsd.edu/chapters/chap11.pdf
[7]
Use the Harvard Referencing style for references in the document.
For more information on when to reference see the following Wikipedia entry.
This page titled 2.7: Solving ODEs with Mathematica- How to find numerical and analytical solutions to ODEs with Mathematica is shared under a CC BY
3.0 license and was authored, remixed, and/or curated by Peter Woolf et al. via source content that was edited to the style and standards of the LibreTexts
platform; a detailed edit history is available upon request.
2.8: FITTING ODE PARAMETERS TO DATA USING EXCEL- USING
REGRESSION TO FIT COMPLEX MODELS IN EXCEL
2.8.1: SECTION 8. FITTING ODE PARAMETERS TO DATA USING EXCEL: USING REGRESSION TO FIT
COMPLEX MODELS IN EXCEL
Note: Video lecture available for this section!
Authors: Anthony Campbell, Victoria Cardine, David Hines, Stephen Kerns /Date Revised: 14th September 2007
Stewards: Jeff Byrd, Khek Ping Chia, John Cruz, Natalie Duchene, Samantha Lyu
2.8.1.1: INTRODUCTION
One of the most important abilities of an engineer is to model measured values of a system. Modeling of a system aids in optimizing
objectives, as well as predicting future operations. This can be especially helpful when changes cannot be immediately implemented. In
general modeling terms, a set of data can always fit a set of variables. However, an excessive number of degrees of freedom usually exist,
creating multiple solutions. Therefore, the physical significance of the variables is an important consideration.
As time proceeds, various values may need to be changed to accommodate variations within the system. The following sections will discuss
the use of Microsoft Excel© as a tool to adjust the values of variables given a set of experimental data. The models presented below take
into account the design equations of the system being modeled, and the rate laws, where applicable.
2.8.1.2: ADDING IN THE SOLVER APPLICATION IN EXCEL 2007

In order to use "Solver" in Excel 2007, one must first add the Excel tool-pack. To add the tool-pack, follow these instructions:
Click on the Microsoft Office symbol at the top left corner of the screen
Click on “Excel options” at the bottom right corner of the drop box
Click on “Add-Ins” in the left hand column
Click on the “Solver Add-In” option at the bottom of the list
Click on “GO” at the bottom right corner
Check the “Solver Add-in” box at the bottom of the list
Click on “OK” at the top of the right column
To access Solver, click on the “Data” tab at the top of the screen. The Solver tool is located in the "Analysis" box on the right side of the
screen.
The user interface of Excel 2007 differs from those of previous versions. The following images illustrate the process just described.
In order to use the Solver function, click on Data and then the Solver button in the Analysis box.
NOTE: If using Microsoft Office 2003, the Solver Application can be added by clicking on the "Tools" tab at the top of the screen. The
same directions can be followed after this step.
NOTE: If you receive an error message while trying to access Solver, it may already be added in Excel. Uninstall Solver (by following the
same steps, but unchecking "Solver Add-in") then reinstall it.
2.8.1.3: ADDING IN THE SOLVER APPLICATION IN EXCEL 2008

Excel 2008 for Macs does not come with the solver add-in. The software can be downloaded at Solver.com
This solver application requires Excel 12.1.2. If this is not the current version of Excel being used, it can be updated by running Microsoft
AutoUpdate or by opening excel, going to help and then clicking on "Check for Updates". Make sure when updating Excel to close all
Microsoft Applications.
This application allows the user to have a solver program that can be opened and used while running Excel. It works the same way as the
solver add-in.
2.8.1.4: USING EXCEL SOLVER TO FIT ODE PARAMETERS
Excel Solver is a tool that can be used to fit function variables to given experimental data. The following process can be used to model data:
Define the data set by entering values into an Excel spreadsheet
Define the model that you want to fit to the data
Define the sum of least squares in one of the spreadsheet blocks
In Solver, have Excel minimize the sum of least squares by varying the parameters in the model
Accurate modeling is dependent on two factors; initial values and verifying results.
Initial values: As stated in the introduction, many data fitting problems can have multiple “solutions”. In numerical methods, given a set of
initial parameters the data will converge to a solution. Different initial values can result in differing solutions. Therefore, if the initial values
are not set by the problem statement, multiple initial guesses should be entered to determine the “best” value for each variable.
Verification of curve fit: When fitting a curve to given data points, it is important to verify that the curve is appropriate. The best way to do
this is to view the data graphically. It is fairly simple to use the Chart tool to graph the data points. Adding a trend line will show the
mathematical relationship between the data points. Excel will generate a function for the trend line and both the function and the R-squared
value can be shown on the plot. This can be done by clicking “Options” next to “Type” of trend and checking the empty boxes that are
labeled “Display equation on chart” and “Display R-squared value on chart”.
There are various types of trendlines that can be chosen for your given data points. Depending on the trend that you may see, you might
want to try a few of the options to obtain the best fit. Sometimes the squared residual value will be better for a certain type of trend, but the
trend is not necessarily “correct”. When you have obtained sample data, be aware of the trend you should be seeing before fitting a certain
trend type to the data points (if this is possible). The polynomial trend type gives you the option to change the order of the equation.
Typically, you will not acquire a squared residual value of zero. This value is a simple analysis of the error between the trend line, and the
actual data. Once the function and squared residual values are generated, you can begin to evaluate the generated solution.
During your solution analysis, it is important to check for variables or parameters that may be unnecessarily influencing the data. Excel can
generate a "best" solution, which is not correct due to background calculations or Excel's misinterpretation of the modeled system. To assess
this possibility, you can graph data generated by using the solved for variables against the given data.
Additional information about the Excel Solver Tool can be found in Excel Modeling.
2.8.1.5: WORKED OUT EXAMPLE 1: MASS BALANCE ON A SURGE TANK

In this example we will walk through the scenario in which you want to model a non-heated surge tank in Microsoft Excel©. This model
would be very critical in a process where smooth levels of process flow are required to maintain product specifications. Therefore, an
engineer would like to know whether they will need to periodically refill or drain the tank to avoid accidents. Let's look at this simple
example:
Problem statement:
Suppose water is being pumped into a 10-ft diameter tank at the rate of 10 ft3/min. If the tank were initially empty, and water were to leave
the tank at a rate dependent on the liquid level according to the relationship Qout = 2h (where the units of the constant 2 are [ft2/min], of Q
are [ft3/min], and of h are [ft]), find the height of liquid in the tank as a function of time. Note that the problem may be, in a way, over
specified. The key is to determine the parameters that are important to fit to the experimental surge tank model.
Solution:
The solution with fitted parameters can be found here: Surge tank solution This Excel file contains this scenario, where given the actual data
(the height of the liquid at different times), the goal is to model the height of the liquid using an ODE. The equations used to model h(t) are
shown and captioned. In order to parameterize the model correctly, Solver was used to minimize the sum of the least square differences at
each time step between the actual and model-predicted heights by changing the parameters' values. The model predicts the actual height
extremely well, with a sum of the least squared differences on the order of 10^-5. A graph showing the modeled and actual height of the
tank as a function of time is included as well.
2.8.1.6: WORKED OUT EXAMPLE 2: FITTING A HEATED SURGE TANK MODEL(A HEATED CSTR) PARAMETERS TO
DATA IN EXCEL
Problem Introduction: In this example we will analyze how we can fit parameters used for modeling a heated CSTR to data stored in
Excel. The purpose of this problem is to learn how to vary ODE parameters in Excel in order to achieve a square difference between the
calculated and experimental data of zero. Some questions to consider are, “How do I model this? Is this a first order ODE? Are there any
coupled ODE's? How many variables need to be considered?" This example will use an Excel file to help you in this process.
Problem Statement: Download the heated CSTR problem file.CSTR Problem Take the experimental data produced from the heated CSTR
model in the Excel spreadsheet*. You are modeling a heated CSTR. It would be wise and is highly recommended to review both the heated
CSTR problem and a wiki on using Excel's solver.
You know this is a heated CSTR...
2. Numerically solve the ODEs through modifying the given wrong set of parameters.
3. While changing the parameters, use Excel to calculate the error between the calculated model and the experimental data.
4. Use the Excel solver function to solve for some of the parameters to minimize error to less than 5%
5. How would your approach be different if our heat transfer area needed to be as small as possible due to capital constraints? Explain.
6. How would your approach be different if our mass feed rate has to be in the range of 3-3.5 kg/s? Explain.
Problem Solution: The solution with fitted parameters can be found here.CSTR Solution
With the model and data given in the example the error must first be defined. Here there error is the squared residual or [(model observation
at time i)-(data observation at time i)]^2 over (all i).
Chemical engineering intuition must be used to reasonably approximate some variables. In real life we would double check that "parameters
normally known, but erroneous" are correct. (Due to problem constraints, we will assume these are accurate, but a true chemical engineer
would check regardless). We must now modify "parameters often needing fitting". For part 3, we can modify the parameters as shown in the
solution.
The original parameters have been modified to produce an error of zero.
If we have restrictions on heat exchange area then we should try to fit the "heat transfer coefficient and area" as low as possible. Likewise, if
we have restrictions on feed rate, we must keep this above 4 while fitting other parameters as best as possible to achieve a small error. These
have certain industrial applications as costing of heat exchangers and limitations on mass flow rates need to be modeled before a reactor
design is implemented.
Side Note: The example above assumes a first order rate law. This is not always the case. To change the equations to fit a new rate law, one
must derive a new dCa/dt equation and dT/dt using their knowledge of reactor engineering. Once the new equations are derived write the
formula in an excel format in our dCa/dt and dT/dt column and then drag it to the end of the column. All other columns should stay the
same.
2.8.1.7: SAGE'S CORNER

A brief narration for a better understanding of Example 1 - Surge Tank Model and tips for effectively picking a model equation can be found
here: video.google.com/googleplayer.swf?docId=70222262469581741
A copy of the slides can be found here:
Unnarrated Slides
An example of fitting an ODE parameter to data using excel - based on a simple radioactive decay rate ODE:
A copy of the slides can be found here: Unnarrated Slides
A copy of the excel file can be found here: Excel File
An example of modelling a reaction rate ODE with rate constant and rate order variable parameters:
video.google.com/googleplayer.swf?docId=-4220475052291185071
A copy of the slides can be found here: Unnarrated Slides
A copy of the excel file can be found here: Excel File
2.8.1.8: REFERENCES
Holbert, K.E. "Radioactive Decay", 2006, ASU Department of Electrical Engineering, Online: October 1, 2007. Available
www.eas.asu.edu/~holbert/eee460/RadioactiveDecay.pdf
This page titled 2.8: Fitting ODE parameters to data using Excel- Using regression to fit complex models in Excel is shared under a CC BY 3.0 license and
2.9: HELPFUL MATHEMATICA SYNTAX- HINTS ON HOW TO USE
MATHEMATICA TO MODEL CHEMICAL PROCESSES
2.9.1: SECTION 9. HELPFUL MATHEMATICA SYNTAX: HINTS ON HOW TO USE MATHEMATICA TO
MODEL CHEMICAL PROCESSES
Title: Useful functions in Mathematica as pertaining to Process controls
2.9.1.1: INTRODUCTION
Mathematica has many different functions that are very useful for controlling chemical engineering processes. The Mathematica syntax of
these functions is often very complicated and confusing. One important resource is the Mathematica Documentation center, which is
accessed through the program help menu. The Documentation Center provides an extensive guide of every function in the program
including syntax and sample code. A "First Five Minutes with Mathematica" tutorial is also present in the menu, which provides a quick
guide to get started with Mathematica. This article will briefly outline several useful Mathematica functions and their syntax.
2.9.1.2: MATHEMATICA BASICS

Here are some basic rules for Mathematica.
Notebook, Cells and Evaluation
The standard file for Mathematica is called a notebook. The notebook consists of two kinds of cells: input and output. The inputs are typed
into one cell and the outputs are displayed in another (See Below).
In order to have Mathematica calculate your inputs, you must evaluate the cells. Do this by pressing Shift+Enter. This will only evaluate
the current cell that your cursor is located. To evaluate all the cells of a notebook press Ctrl+A, to select all, and then press Shift+Enter.
Parentheses vs. Brackets
Brackets, [], are used for all functions in Mathematica. The only thing parentheses, (), are used for is to indicate the order of operations.
Equality, = vs. ==
The single equal sign (=), is used to assign a value or a function to a variable. The double equal sign (==) is used to set two values equal to
each other, such as solving for a value, or to test equalities. When testing equalities, Mathematica will output 'True' or 'False.'
If you recieve 'True' or 'False' as an output when you are expecting it, you likely used a == somewhere that you should have used a =.
Semicolon Use
Placing a semicolon after a line of code will evaluate the expression; however, there will be no output shown.
Mathematica Functions
All Mathematica functions start with a capital letter, eg. Sin[x], Exp[x], Pi or Infinity.
Assigning and Inserting Variables or Parameters
To assign values to variables, use -> rather than an equals sign. For example,
To insert a set of parameters to a function use /. This symbol applies a rule or list of rules in an attempt to transform each subpart of an
expression.
For example, if you want to enter the above parameters into the expression y = Kc * x + Ca0, enter the following in Mathematica:
Variables are case sensitive. For example, X is different than x. Also, if you define a variable, it will remain defined until you redefine it as
something else or clear it. You can clear a variable, for example, x, by using the Clear[x] command. You can also quit the Kernel, which will
clear all variables.
Forcing a Numerical Expression
To force Mathematica into giving a numerical solution, use the function N[expression]. For example,
Another method to do this is to place a decimal place in any of your numbers (i.e. "5." instead of just "5")
2.9.1.3: INTEGRATION
Mathematica can do some very complex integrations with fairly little coding. Integration can be very helpful for many engineering
applications, but when the functions become very complex they become increasingly difficult to integrate. Mathematica will integrate any
function that is possible to integrate.
To integrate a function f(x) the Mathematica, the command is Integrate[].
For an indefinite integral the correct coding is Integrate[f(x),x] where x is the independent variable.
For a definite integral with limits of y and z the correct coding is Integrate[f(x),{x,y,z}]
For example:
We can integrate a function like the one below:
f(x) = Sin(5 * x / 8) * x
We can find the indefinite integral as seen here:
We can find the definite integral from 0 to π as seen here:
2.9.1.4: SOLVER
Mathematica's Solve function attempts to solve an equation or set of equations for the specified variables. The general syntax is Solve[eqns,
vars]
Equations are given in the form of left hand side == right hand side.
A single variable or a list of variables can be also be specified. When there are several variables, the solution is given in terms of a list of
rules. For example,
When there are several solutions, Solve gives a list of them. For example,
Solve will not always be able to get explicit solutions to equations. It will give the explicit solution if it can, then give a symbolic
representation of the remaining solutions in terms of root objects. If there are few symbolic parameters, you can the use NSolve to get
numerical approximations to the solutions.
2.9.1.5: PLOTTING
Mathematica contains a very powerful plotting interface and this overview will go over the basics of plotting. Using Mathematica's Plot[ ]
function there is a basic syntax to be used:
Plot[function, {variable, variable_min, variable_max}]
For example, the command "Plot[y=x, {x,0,10}]" will output a graph of the function y=x for a range of x from 0 to 10.
The function could also be previously defined and then referenced in the Plot [ ] command, this functionality provides a lot of simplicity
when referencing multiple functions and helps prevent a typo from affecting your Plots.
y[x]==x
Plot[y[x],{x,0,10}]
To plot multiple functions on the same plot, use an array of functions using curly brackets { }, in place of the single one:
Plot[{y=x,y=2*x}, {x,0,10}]
Another useful plotting option in Mathematica is the ListPlot function. To use this function, the following formatting should be used:
ListPlot[Table[f(x),{x,min,max}]
This is obviously only the basics, Mathematica provides many more options including colors, legends, line styles, etc. All of these extra
features are well documented at Wolfram's Reference Site.
An additional piece of information: Type the following command on the line after the Plot command to determine the max value plotted. It
will display the number associated with the absolute max of the function within the plotting range.
Max[Last/@Level[Cases[%,_Line,Infinity],{-2}]]
2.9.1.6: SOLVING ODES

Within Mathematica, you have the option of solving ODEs numerically or by having Mathematica plot the solution.
Explicit solution:
The built-in function DSolve solves a differential equation or a list of differential equations for either one indepedent variable or a list of
independent variables, which solves a partial differential equation. The most basic syntax for the function DSolve is DSolve[eqn,y,x]. For
example,
The syntax for solving a list of differential equations is DSolve[{eqn1,eqn2,...},{y1,y2,...},x], while the syntax for solving a partial
differential equation is DSolve[{eqn, y, {x1,x2,...}]. For example,
Note that differential equations must be stated in terms of derivatives, such as y'[x]. Boundary conditions can also be specified by giving
them as equations, such as y'[0]==b.
Plotted Solution:
Solutions given by DSolve sometimes include integrals that cannot be carried out explicitly, thus NDSolve should be used. NDSolve finds a
numerical solution to an ordinary differential equation within a specified range of the independent variable. The result will be given in terms
of InterpolatingFunction objects, and thus the solution should be plotted. For example,
Note that a place holder was defined to store the NDSolve solution in order for the solution to be plotted. The general syntax for plotting is
Plot[Evaluate[y[x] /. s], {x, xmin, xmax}], where s is said placeholder. The plot below shows the solution to this example.
2.9.1.7: ANOTHER EXAMPLE SOLVING ODES

In this example the water level of a tank is controlled by adjusting the feed rate of water into the tank. The system is modeled by the
following equation:
There is a PID controller on the valve but is not important in explaining the Mathematica syntax. Consider it just a different way to define
Fin.
Parameters are assumed to be Vset = 10, Kc = 10, tauI = 0.1, and tauD = 1
For this problem the differential equation needs to be modeled with the given parameters in Mathematica. The problem asks what the steady
state of the system is. This is solved for by substituting in a PID expression for Fin (not important in this context), setting the derivatives
equal to zero, and solving for V and Xi(a made up variable for the PID controller). The following Mathematica code does all said things.
The final line shows that the steady state values for the problem are V = Vset or 10 units of volume, and Xi (the made up variable for the
PID controller) has a steady state value of 3.9
In this example there is syntax for defining parameters, applying those parameters to equations, setting up equations, and solving equations
for steady state values. There is also syntax for plotting multiple equations on the same graph and how to evaluate multiple functions over
time periods. Overall very useful for solving Controls Engineering problems in Mathematica.
2.9.1.8: MATRIX
It is very simple to create a matrix in Mathematica. Simply pick a name for the matrix (for this example we will use a) then simply enter
each row inside of curly brackets {} with the individual values separated by commas, and the the row separated by a comma. The function
MatrixForm[] will display a matrix in its matrix form rather than as a list of numbers.
EX: a =
ParseError: EOF expected (click for details)
MatrixForm[a] displays:
523
7 9 11
8 13 9
You can refer to any specific position in the matrix by using the system:
matrix[[row],[column]]
EX: For the matrix above a[[1],[2]] a value of 2 would be returned.
2.9.1.9: JACOBIANS AND EIGENVALUES

The Jacobian is a helpful intermediate tool used to linearize systems of ODEs. Eigenvalues can be used with the Jacobian to evaluate the
stability of the system. Mathematica can be used to compute the Jacobian and Eigenvalues of a system by the following steps.
1. Enter the system as individual equations
entered as
eq1 = 3*x^2+2y
eq2 = x^2-3*y^2
2. Compute the Jacobian with Mathematica's Derivate function
computed by
Jac =
ParseError: "]" expected (click for details)
3. Use Mathematica's Eigenvalues function to evaluate the stability of the system

Begin by entering a desired value of {x,y} into the matrix (i.e. steady state values of x and y).
Jac =
Then call the Eigenvalue function.
Eigenvalues[Jac]
Result: λ={-1.82, 0.85}
The Jacobian can be computed as long as the matrix is square. Many other function types are acceptable when entering eq1 and eq2,
including trigonometric and exponential functions.
2.9.1.10: USER-DEFINED FUNCTIONS

Mathematica also supports user-defined functions. Defining your own functions helps to simplify code, especially when lengthy expressions
are used consistently throughout the code. The generic Mathematica syntax for defining a function is as follows:
f[a_,b_,c_]:=(a+b)/c
Here, the function name is "f" with three arguments: a, b, and c. The syntax ":=" assigns the expression on the right hand side to the function
"f." Make sure that each of the arguments in the brackets on the left hand side are followed by underscores; this identifies the variables
being called by the function as local variables (ones that are only used by the function itself). For this example, the function "f" will take
three arguments (a, b, and c) and calculate the sum of a and b, and then divide by c.
After defining the function, you can use the function throughout your code using either real numbers, as the first example shows, or any
variables. This is done by calling the function with the correct number of arguments. When referencing variables, if the variable already has
an assigned value the value will be used, if it does not, Mathematica will display your function's results with the variables displayed. A few
examples are shown below.
2.9.1.11: NONLINEAR REGRESSION

Mathematica can be used to solve multiple variable, nonlinear regression problems. This form of regression can be useful to create models
of how temperature or pressure sensors depend on different valves in a system. The NonLinearRegress function in Mathematica can be used
to find a best-fit solution to a user defined functional form.
The NonLinearRegress function is not part of the standard package in Mathematica. For this reason, this function must be loaded into
Mathematica before an attempt to use.
After this step, the NonlinearRegress is loaded and can be used to solve for a set of data. To use the this function, enter
NonlinearRegress[data, expression, parameters, variables]. The best fit values for the parameters will be returned along with statistical
information about how well this form fits the given data. An example of Nonlinear Regression in Mathematica is shown in the image below.
2.9.1.12: PROBABILITY DENSITY FUNCTION
Often in controls, it is necessary to calculate the probability that a certain event will occur, given a particular known distribution. In these
scenarios, it is helpful to know how to use the built-in probability density function in Mathematica. The syntax for using the probability
density function to calculate the probability of observing a value x for a generic distribution dist is shown below.
PDF[dist,x]
A variety of distributions can be used in conjunction with the probability density function. These distributions include the normal
distribution and the binomial distribution.
Normal Distribution
For a normal distribution with a mean μ and a standard deviation σ, use the following syntax to specify the distribution.
NormalDistribution[μ,σ]
To calculate the probability of observing a value x given this distribution, use the probability density function as described previously.
PDF[NormalDistribution[μ,σ],x]
To calculate the probability of observing a range of values from x1 to x2 given this distribution, use the following syntax.
NIntegrate[PDF[NormalDistribution[μ,σ],x],{x,x1,x2]
For the theory behind the probability density function for a normal distribution, click here.
Binomial Distribution
For a binomial distribution with a number of trials n and a success probability p, use the following syntax to specify the distribution.
BinomialDistribution[n,p]
To calculate the probability of observing a successful trial exactly k times given this distribution, use the probability density function as
described previously.
PDF[BinomialDistribution[n,p],k]
To calculate the probability of observing between a range of k1 number of successes and k2 number of successes given this distribution, use
the following syntax.
NIntegrate[PDF[BinomialDistribution[n,p],k],{k,k1,k2]
Alternatively, binomial distribution can be found using a user defined function. Define the equations as follow: binom[nn_,kk_,pp_]:=
nn!/(kk! (nn-kk)!) pp^kk (1-p)^(nn-kk)
where:
number of independent samples = nn
number of events = kk
probability of the event = pp
To find the odds of getting 5 heads out of 10 coin tosses,assuming the probability of head = 0.5, substitute the values for nn,kk and pp as
following: binom[10,5,0.5]
Using the same scenario, but instead for the odds of getting 5 or more heads in the 10 tosses, the function can be used as:
Sum[binom[10,i,0.5],{i,5,10}]
For the theory behind the probability density function for a binomial distribution, click here.
2.9.1.13: MANIPULATE FUNCTION

Often in controls, it is useful to manipulate certain variables dynamically and see the result in a plot or table. This is particularly useful in
setting the constants in tuning PID controllers. The Manipulate function in Mathematica lets you do this. You can use the manipulate
function in conjunction with the Plot command. The general syntax for using the manipulate function is as follows:
Manipulate[expr,{u,umin,umax}]
Where expr is the function that you want to add controls to and u is the variable you wish to manipulate.
Manipulate a function:
You can change the variable x in the function s. Moving the bar lets you change the variable x.
Manipulate a plot:
You can also manipulate a plot. The slider bar lets change the variable n and the plot changes accordingly.
The Mathematica file used to create this example, can be downloaded below.
Manipulate Command Example
For more information regarding manipulate, please see the Wolfram's Reference Site.
[Note - need to add manipulate for the Bode plots - R.Z.]
2.9.1.14: ADDING GRIDLINES
2.9.1.15: PLOTTING LOG-LOG GRAPHS
2.9.1.16: OTHER USEFUL MATHEMATICA TIPS
[It would be good to replace this images with actual text - R.Z.]
This page titled 2.9: Helpful Mathematica Syntax- Hints on how to use Mathematica to model chemical processes is shared under a CC BY 3.0 license and

Chapter Overview: 2: Modeling Basics

Uploaded by

Copyright:

Available Formats

Chapter Overview: 2: Modeling Basics

Uploaded by

Document Information

Original Title

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

Chapter Overview: 2: Modeling Basics

Uploaded by

Copyright:

Available Formats

CHAPTER OVERVIEW

2.1.2: PREREQUISITE INFORMATION REGARDING A PROCESS

2.1.3: STEP-BY-STEP METHOD FOR DESCRIBING CONTROLS AND THEIR PURPOSE

2.1.4: ALTERNATIVE METHOD OF VERBAL MODELING

2.1.5: THE BARKEL METHOD OF VERBAL MODELING

2.1.6: COMMON ERRORS

2.1.7: WORKED OUT EXAMPLE 1

2.1.9: SAGE'S CORNER

2.2.3: THE PROCESS

2.2.5: EXAMPLE 1: BLENDER

Figure 2: Blender with controls

Monotonic Increase : If x≥y, then f(x)≥f(y)

Monotonic Decrease : If x≥y, then f(x)≤f(y)

2.3.2: INCIDENCE GRAPHS

This incidence graph indicates that:

This incidence graph indicates that:

This incidence graph indicates that:

2.3.3: CONSISTENT, INCONSISTENT, AND PARTIALLY CONSISTENT

2.3.3.1: CONSISTENT GRAPHS

Paths(s) Sign series

All pathways leading from node 1 to node 3 -

Paths(s) Sign series

2.3.3.2: INCONSISTENT GRAPHS

Paths(s) Sign series

Paths(s) Sign series

All pathways leading from node 1 to node 2 -

2.3.3.3: PARTIALLY CONSISTENT GRAPHS

Paths(s) Sign series

All pathways leading from node 1 to node 4 -

2.3.4: WORKED OUT EXAMPLE 1

2.3.5: WORKED OUT EXAMPLE 2

2.3.6: WORKED OUT EXAMPLE 3

2.3.8: CONTRIBUTORS AND ATTRIBUTIONS

2.4.1.1: LOGICAL METHODS

3) Click solve to find the solution.

2.4.1.3: NONLINEAR REGRESSION

2.4.1.3.1: ANOVA ANALYSIS

2.4.1.4: RANDOM NUMBER SAMPLING

Microsoft Excel 2007:

2.4.1.5: WORKED OUT EXAMPLE 1

3. Click the Min radio button

2.4.1.6: WORKED OUT EXAMPLE 2

2.4.1.7: ACCESSING EXCEL TOOLS ON WINDOWS

2.4.2: CONTRIBUTORS AND ATTRIBUTIONS

2.5.2: GENERAL TYPES OF NOISE

Microsoft Excel Formulas for 3 Types of Noise

2.5.3: COLORS OF NOISE

2.5.3.1: WHITE NOISE

image from en.Wikipedia.org/wiki/White_noise on 9/18/2006

2.5.3.2: PINK NOISE

image from en.Wikipedia.org/wiki/Pink_noise on 9/18/2006

2.5.3.3: BROWN NOISE

image from en.Wikipedia.org/wiki/Brown_noise on 9/18/2006

2.5.3.4: BLUE NOISE

image from en.Wikipedia.org/wiki/Blue_noise#Blue_noise on 12/11/2008

2.5.3.5: PURPLE NOISE

2.5.3.6: MODELING COLORS OF NOISE

P S D = √ (∑ x(t) cos(wt)) + (∑ x(t) sin(wt))

2.5.4: POPS, SNAPS & CRACKLES