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
appl = pd.read_csv('http://apmonitor.com/che263/uploads/Main/AAPL.csv')
goog = pd.read_csv('http://apmonitor.com/che263/uploads/Main/GOOG.csv')
#xom = pd.read_csv('http://apmonitor.com/che263/uploads/Main/XOM.csv')
# create dictionary of stocks
s = dict([('Apple',appl),('Google',goog)]) #,('ExxonMobil',xom)])

# print column headers and starting rows (5 is default)
print('Apple Data')
print(s['Apple'].head())

# print column headers and ending close price (4 rows)
print('Google Data')
print(s['Google']['Close'].tail(4))

# 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()
sty = dict([('Apple','r--'),('Google','b:'),('ExxonMobil','k-')])
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
x = pd.read_csv('http://apmonitor.com/che263/uploads/Main/tclab.txt')

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

# 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)
url = 'http://apmonitor.com/che263/uploads/Main/heart_rate.txt'
x = pd.read_csv(url)

# print first rows
print('Data')
print(x.head())

# 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)
url = 'http://apmonitor.com/che263/uploads/Main/heart_rate.txt'
x = pd.read_csv(url)

# print first rows
print('Data')
print(x.head())

# 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
url = 'http://apmonitor.com/che263/uploads/Main/dynamics.txt'
x = pd.read_csv(url)

# print first rows
print('Data')
print(x.head())

# 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
def dCadt(t,Ca):
    return -k * Ca
Ca_m4 = odeint(dCadt,t,Ca0)

# 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')

Python Assignments

MathCAD Assignments