Main

## Programming Assignments

#### Excel Assignments

Complete the following assignments with Microsoft Excel, Google Sheets, or another similar spreadsheet program. Microsoft Excel solutions are shown for each problem but should only be used as a learning resource, not to merely complete the assignment. Python solutions are shown as examples of programming but are not required.

• Assignment #1
• Conditionals
• Functions
• Excel Solution 2, Solution 3
# constants
L      = 2       # m
Tplate = 343     # K
v      = 1.45    # m/s
Twater = 294     # K
mu     = 9.79e-4 # Pa*s
rho    = 998     # kg/m^3
k      = 0.601   # W/m-K
cp     = 4.18e3  # J/kg-K

# derived quantities
Re = rho*L*v/mu
Pr = mu*cp/k
Nu = 0.332 * Pr**(1.0/3.0) * Re**(1.0/2.0)
h = Nu * k / L
q = h * (Tplate-Twater)

print('Rate of Heat Transfer')
print(str(q)+' W/m^2')
import numpy as np
import matplotlib.pyplot as plt

# calculate initial moles of N2
T = 298.15  # K
P = 1.25    # atm
V = 4000    # L
Rg = 0.0821 # L*atm/mol/K
n_N2 = P*V/(Rg*T) # moles of N2
n_N2 = n_N2 * 0.25  # moles with 25% remaining

# calculate N2 pressure at different V and T
n = 20 # grid points
V = np.linspace(500,1000,n)
T = np.linspace(100,600,n)
# create meshgrid
Vm,Tm = np.meshgrid(V,T)
# calculate pressure
Pm = n_N2 * Rg * Tm / Vm

# plot results
plt.figure()
plt.contourf(Vm,Tm,Pm,cmap='RdBu_r')
plt.colorbar()
CS2 = plt.contour(Vm,Tm,Pm,[1.0,2.5],colors='k')
plt.clabel(CS2, inline=1, fontsize=10)
plt.xlabel('Volume (L)')
plt.ylabel('Temperature (K)')
plt.show()
import pandas as pd
import matplotlib.pyplot as plt
import numpy as np

# import April 2018 data or get new data from finance.yahoo.com
# create dictionary of stocks

# print column headers and starting rows (5 is default)
print('Apple Data')

# print column headers and ending close price (4 rows)

# basic data statistics
for i in s:
print('Stock: ' + i)
print(' max   : ' + str(max(s[i]['Close'])))
print(' min   : ' + str(min(s[i]['Close'])))
print(' stdev : ' + str(np.std(s[i]['Close'])))
print(' avg   : ' + str(np.mean(s[i]['Close'])))
print(' median: ' + str(np.median(s[i]['Close'])))

# plot data
plt.figure()
ni = 0
for i in s:
mc = max(s[i]['Close'])
plt.plot(s[i]['Date'],s[i]['Close']/mc,sty[i],linewidth=3,label=i)
plt.plot(s[i]['Date'],s[i]['High']/mc,sty[i],linewidth=1)
plt.plot(s[i]['Date'],s[i]['Low']/mc,sty[i],linewidth=1)
plt.xticks(rotation=90)
plt.legend(loc='best')
plt.show()
import pandas as pd
import matplotlib.pyplot as plt

# import data
# Time (sec), Heater 1, Heater 2, Temperature 1, Temperature 2

# print column headers and starting 10 rows (5 is default)
print('Data')

# plot data
plt.figure()
plt.subplot(2,1,1)
plt.title('Temperature Control Lab')
plt.plot(x['Time (sec)'],x['Temperature 1'],'r--')
plt.plot(x['Time (sec)'],x['Temperature 2'],'b-')
plt.ylabel('Temp (degC)')
plt.legend(loc='best')

plt.subplot(2,1,2)
plt.plot(x['Time (sec)'],x['Heater 1'],'r--')
plt.plot(x['Time (sec)'],x['Heater 2'],'b-')
plt.ylabel('Heater (%)')
plt.xlabel('Time (sec)')
plt.legend(loc='best')
plt.show()
import math as m

# part a
help(m.cos)
y = m.cos(0.5)
print('cos(0.5): ' + str(y))

# part b
y = m.sin(30.0*(m.pi/180.0))
print('sin(30 deg): ' + str(y))

# part c
y = m.tan(m.pi/2.0)
print('tan(pi/x): ' + str(y))

# part d
x = 5.0
y = max(2.0*m.sqrt(x),(x**2)/2.0,\
(x**3)/3.0,(x**2+x**3)/5.0)
print('max value: ' + str(y))

# part e
x = 25
y = m.factorial(x)
print('25!: ' + str(y))

# part f
x = 0.5
# if..else statement
if x<1.0:
y = x**2
else:
y = m.sin(m.pi*x/2.0)
# same but one line
y = x**2 if x<1.0 else m.sin(m.pi*x/2.0)
print('if statement result: ' + str(y))

# part g
x = 4.999
y = m.floor(x)
print('floor(4.999): ' + str(y))
# method #1: NumPy
from numpy.linalg import solve
A = [[11.0, 3.0, 0.0,  1.0,  2.0],\
[ 0.0, 4.0, 2.0,  0.0,  1.0],\
[ 3.0, 2.0, 7.0,  1.0,  0.0],\
[ 4.0, 0.0, 4.0, 10.0,  1.0],\
[ 2.0, 5.0, 1.0,  3.0, 14.0]]
b = [45.0, 30.0, 15.0, 20.0, 92.0]
x = solve(A,b)
print('NumPy Solution')
print(x)

# method #2: Gekko
from gekko import GEKKO
m = GEKKO()
x1,x2,x3,x4,x5 = [m.Var() for i in range(5)]
m.Equation(11*x1+3*x2+x4+2*x5==45)
m.Equation(4*x2+2*x3+x5==30)
m.Equation(3*x1+2*x2+7*x3+x4==15)
m.Equation(4*x1+4*x3+10*x4+x5==20)
m.Equation(2*x1+5*x2+x3+3*x4+14*x5==92)
m.solve(disp=False)
print('Gekko Solution')
print(x1.value)
print(x2.value)
print(x3.value)
print(x4.value)
print(x5.value)
# method #1: NumPy
from scipy.optimize import fsolve
def f(z):
x,y=z
f1 = 2*x**2+y**2-1
f2 = (0.5*x-0.5)**2+2.0*(y-0.25)**2-1
return [f1,f2]
x,y = fsolve(f,[1,1])
print('NumPy Solution')
print(x,y)

# method #2: Gekko
from gekko import GEKKO
m = GEKKO()
x,y = [m.Var(value=1) for i in range(2)]
m.Equation(2*x**2+y**2==1)
m.Equation((0.5*x-0.5)**2+2.0*(y-0.25)**2==1)
m.solve(disp=False)
print('Gekko Solution')
print(x.value)
print(y.value)
import numpy as np
from scipy.optimize import fsolve

# constants
TC = 77 # degC
P = 1.0 # bar
a = 2.877e8 # cm^6 bar K^0.5 / mol^2
b = 60.211  # cm^3 / mol
Rg = 83.144598 # cm^3 bar / K-mol
# derived quantities
TK = TC+273.15 # K

# method #1: NumPy
def f(V):
return P - Rg*TK/(V-b)+a/(np.sqrt(TK)*V*(V+b))
V_liq = fsolve(f,82)    # Liquid root
V_vap = fsolve(f,28600) # Vapor root
print('NumPy Solution')
print(V_liq,V_vap)

# method #2: Gekko
from gekko import GEKKO
m = GEKKO()
V = m.Var(value=[82,28600])
m.Equation(P==Rg*TK/(V-b)-a/(m.sqrt(TK)*V*(V+b)))
m.options.IMODE=2
m.solve(disp=False)
print('Gekko Solution')
print(V.value)
from gekko import GEKKO

# constants
y1 = 0.33
y2 = 1.0-y1
P = 120.0 # kPa

# Antoine constants
# Benzene
ac1 = [13.7819, 2726.81, 217.572]
# Toluene
ac2 = [13.9320, 3056.96, 217.625]

# gekko model
m = GEKKO()
T = m.Var(value=100)
x1,x2 = [m.Var(value=0.5,lb=0,ub=1) for i in range(2)]
# vapor pressure
Psat1 = m.Intermediate(m.exp(ac1[0]-ac1[1]/(T+ac1[2])))
Psat2 = m.Intermediate(m.exp(ac2[0]-ac2[1]/(T+ac2[2])))
# Raoult's law
m.Equation(y1*P==x1*Psat1)
m.Equation(y2*P==x2*Psat2)
m.Equation(x1+x2==1)
m.options.IMODE=1
m.solve(disp=False)
print('Gekko Solution')
print('x1: ' + str(x1.value))
print('x2: ' + str(x2.value))
print('T:  ' + str(T.value))
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit
from sklearn.metrics import r2_score

# import data
# Time (sec),Heart Rate (BPM)

# print first rows
print('Data')

# extract vectors
t = x['Time (sec)'].values
ym = x['Heart Rate (BPM)'].values

# define function for fitting
def bpm(t,c0,c1,c2,c3):
return c0+c1*t-c2*np.exp(-c3*t)

# find optimal parameters
p0 = [100,0.01,100,0.01]        # initial guesses
c,cov = curve_fit(bpm,t,ym,p0)  # fit model

# print parameters
print('Optimal parameters')
print(c)

# calculate prediction
yp = bpm(t,c[0],c[1],c[2],c[3])

# calculate r^2
print('R^2: ' + str(r2_score(ym,yp)))

# plot data and prediction
plt.figure()
plt.title('Heart Rate Regression')
plt.plot(t/60.0,ym,'r--',label='Measured')
plt.plot(t/60.0,yp,'b-',label='Predicted')
plt.ylabel('Rate (BPM)')
plt.xlabel('Time (min)')
plt.legend(loc='best')
plt.show()
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from gekko import GEKKO
from sklearn.metrics import r2_score

# import data
# Time (sec),Heart Rate (BPM)

# print first rows
print('Data')

# extract vectors
t = np.array(x['Time (sec)'])
ym = np.array(x['Heart Rate (BPM)'])

# GEKKO model
m = GEKKO()

# parameters
tm = m.Param(value=t)
c0 = m.FV(value=100)
c1 = m.FV(value=0.01)
c2 = m.FV(value=100)
c3 = m.FV(value=0.01)
c0.STATUS=1
c1.STATUS=1
c2.STATUS=1
c3.STATUS=1

# variables
bpm = m.CV(value=ym)
bpm.FSTATUS=1

# regression equation
m.Equation(bpm==c0+c1*tm-c2*m.exp(-c3*tm))

# regression mode
m.options.IMODE = 2

# optimize
m.solve()

# print parameters
print('Optimal parameters')
print(c0.value[0])
print(c1.value[0])
print(c2.value[0])
print(c3.value[0])

# calculate r^2
print('R^2: ' + str(r2_score(ym,bpm)))

# plot data and prediction
plt.figure()
plt.title('Heart Rate Regression')
plt.plot(t/60.0,ym,'r--',label='Measured')
plt.plot(t/60.0,bpm,'b-',label='Predicted')
plt.ylabel('Rate (BPM)')
plt.xlabel('Time (min)')
plt.legend(loc='best')
plt.show()
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit
from sklearn.metrics import r2_score

# import data
# time (min),y

# print first rows
print('Data')

# extract vectors
t = x['time (min)'].values
ym = x['y'].values

# define function for fitting
def yfcn(t,tau,theta):
n = len(t)
res = np.zeros(n)
for i in range(n):
if i>=theta:
res[i] = 5.0*(1.0-np.exp(-(t[i]-theta)/tau))
return res

# find optimal parameters
c,cov = curve_fit(yfcn,t,ym)

# print parameters
print('Optimal parameters')
tau = c[0]
theta = c[1]
print(c)

# calculate prediction
yp = yfcn(t,tau,theta)

# calculate r^2
print('R^2: ', r2_score(yp,ym))

# plot data and prediciton
plt.figure()
plt.title('Step Test Regression')
plt.plot(t,ym,'ro',label='Measured')
plt.plot(t,yp,'b-',label='Predicted')
plt.ylabel('Response')
plt.xlabel('Time (min)')
plt.legend(loc='best')
plt.show()
import numpy as np
import matplotlib.pyplot as plt

# generate 1000 random numbers
#  with Poisson distribution and lambda=1
n = 1000
lam = 1
x = np.random.poisson(lam,n)

# count number in each bin
bins=[0,1,2,3,4,5,6]
hist, _ = np.histogram(x, bins)

# plot histogram data
plt.bar(bins[0:-1],hist,label='1000 samples')
plt.xlabel('bin')
plt.ylabel('count')
plt.title('Poisson Distribution (lambda=1)')
plt.legend(loc='best')
plt.show()
import numpy as np
import matplotlib.pyplot as plt
from gekko import GEKKO
from scipy.integrate import odeint

# number of time points
n = 15
# final time
tf = 7.0
# initial concentration
Ca0 = 5.0
# constants
k = 1 # 1/s

# method #1 Analytical solution
# Ca(t) = Ca(0) * exp(-k*t)
t = np.linspace(0,tf,n)
Ca_m1 = Ca0 * np.exp(-k*t)

# method #2 Euler's method
Ca_m2 = np.empty(n)
Ca_m2[0] = Ca0 # kmol/m^3
for i in range(1,n):
dt = t[i] - t[i-1]
Ca_m2[i] = Ca_m2[i-1] - k * Ca_m2[i-1] * dt

# method #3: GEKKO solution
# create new gekko model
m = GEKKO()
# integration time points
m.time = t
# variables
Ca = m.Var(value=Ca0)
# differential equation
m.Equation(Ca.dt()==-k*Ca)
# set options
m.options.IMODE = 4 # dynamic simulation
m.options.NODES = 3 # collocation nodes
# simulate ODE
m.solve()

# method #4: ODEINT from SciPy
return -k * Ca

# plot results
plt.figure(1)
plt.plot(t,Ca_m1,'r-',label='Ca (Analytical)')
plt.plot(t,Ca_m2,'ko',label='Ca (Euler)')
plt.plot(t,Ca,'b--',label='Ca (GEKKO)')
plt.plot(t,Ca,'ys',label='Ca (ODEINT)')
plt.xlabel('Time (sec)')
plt.ylabel(r'$C_a (kmol/m^3)$')
plt.legend(loc='best')
plt.show()
import numpy as np
import matplotlib.pyplot as plt
from gekko import GEKKO
from scipy.integrate import odeint

# number of time points
n = 15
# final time
tf = 7.0
# initial concentration
Ca0 = 5.0
# constants
k = 1 # 1/s

# method #1 Analytical solution
# Ca(t) = Ca(0) * exp(-k*t)
t = np.linspace(0,tf,n)
Ca_m1 = Ca0 * np.exp(-k*t)

# method #2 Euler's method
t2 = np.arange(0,tf,0.5)
n = len(t2)
Ca_m2 = np.empty_like(t2)
Ca_m2[0] = Ca0 # kmol/m^3
for i in range(1,n):
dt = t2[i] - t2[i-1]
Ca_m2[i] = Ca_m2[i-1] - k * Ca_m2[i-1] * dt

t3 = np.arange(0,tf,1.5)
n = len(t3)
Ca_m3 = np.empty_like(t3)
Ca_m3[0] = Ca0 # kmol/m^3
for i in range(1,n):
dt = t3[i] - t3[i-1]
Ca_m3[i] = Ca_m3[i-1] - k * Ca_m3[i-1] * dt

t4 = np.arange(0,tf,2.1)
n = len(t4)
Ca_m4 = np.empty_like(t4)
Ca_m4[0] = Ca0 # kmol/m^3
for i in range(1,n):
dt = t4[i] - t4[i-1]
Ca_m4[i] = Ca_m4[i-1] - k * Ca_m4[i-1] * dt

# plot results
plt.figure(1)
plt.plot(t,Ca_m1,'r-',label='Ca (Analytical)')
plt.plot(t2,Ca_m2,'k.-',label='Ca (Euler 0.5)')
plt.plot(t3,Ca_m3,'bo-',label='Ca (Euler 1.5)')
plt.plot(t4,Ca_m4,'y--',label='Ca (Euler 2.1)')
plt.xlabel('Time (sec)')
plt.ylabel(r'$C_a (kmol/m^3)$')
plt.legend(loc='best')
plt.show()
import numpy as np
import matplotlib.pyplot as plt
from gekko import GEKKO

# final time
tf = 3.0

# constants
k1 = 1.0 # L/mol-s
k2 = 1.5 # L/mol-s

# GEKKO solution
# create new gekko model
m = GEKKO()
# integration time points
m.time = np.arange(0,tf+0.01,0.2)
# variables
Ca = m.Var(value=1.0)
Cb = m.Var(value=1.0)
Cc = m.Var(value=0.0)
Cd = m.Var(value=0.0)
S = m.Var(value=1.0)
# differential equations
m.Equation(Ca.dt()==-k1*Ca*Cb)
m.Equation(Cb.dt()==-k1*Ca*Cb-k2*Cb*Cc)
m.Equation(Cc.dt()== k1*Ca*Cb-k2*Cb*Cc)
m.Equation(Cd.dt()== k2*Cb*Cc)
m.Equation(S==Cc/(Cc+Cd))
# set options
m.options.IMODE = 4 # dynamic simulation
m.options.NODES = 3 # collocation nodes
# simulate ODE
m.solve()

# plot results
plt.figure(1)
plt.subplot(2,1,1)
plt.plot(m.time,Ca,'r-',label='Ca',linewidth=2.0)
plt.plot(m.time,Cb,'k.-',label='Cb',linewidth=2.0)
plt.plot(m.time,Cc,'b--',label='Cc',linewidth=2.0)
plt.plot(m.time,Cd,'y:',label='Cd',linewidth=3.0)
plt.ylabel('Conc (mol/L)')
plt.legend(loc='best')

plt.subplot(2,1,2)
plt.plot(m.time,S,'k-',label='Selectivity',linewidth=2.0)
plt.xlabel('Time (sec)')
plt.legend(loc='best')
plt.show()
• Assignment #6
• VBA Macros
• Excel Solution 1, Solution 2
import numpy as np

r = 5     # m
h = 10    # m
F = 15    # m^3/min
t = 180   # min

V_tank = np.pi * r**2 * h  # m^3
V_crude_oil = F * t # m^3

if V_crude_oil > V_tank:
print('Tank Overfilled by ' \
+ str(V_crude_oil-V_tank) + ' m^3')
else:
print('Not Overfilled')