This preview has intentionally blurred sections. Sign up to view the full version.
View Full DocumentThis preview has intentionally blurred sections. Sign up to view the full version.
View Full DocumentThis preview has intentionally blurred sections. Sign up to view the full version.
View Full DocumentThis preview has intentionally blurred sections. Sign up to view the full version.
View Full DocumentThis preview has intentionally blurred sections. Sign up to view the full version.
View Full Document
Unformatted text preview: Norse <63, Kinks The University of Queensland
School of Engineering MECH2700 Engineering Analysis I (2005) Exercise Sheet 8 1. A hot—wire anemometer calibration gave the following voltagevelocity pairs. U(m/s) 91.08 85.04 78.54 68.51 58.35 49.63 41.45 34.32
E(m/s) 0.8560 0.8414 0.8375 0.8160 0.7977 0.7762 0.7535 0.7347 The model function E = a0 + a1U + 02W is to be ﬁtted to the calibration data
using the “least—squares” error critereon. 0 Write the design matrix for the given model equation and data points. 0 Formulate the normal equations in matrix form and solve to obtain values for
' the coeﬂicients a0, a1 and a2. 0 Plot E versus U over the range 35 m/s 3 U S 90 m/s for both the original
data and the model function. 2. Solve x = W for N = 0.6 by minimizing the function f (at) = [N — m3]2. Use the “Golden Section” search method, implemented as a Python function. 3. The sampled data for a parameter‘ﬁtting exercise (6.9.1) in the Schilling & Harris
text is given in the following table. It speciﬁes the measured temperature in the flow downstream of a heat exchanger that has a step change in its input at t = 0. 15;, (min) yk 0C 15;, (min) yk °C 0 —0.036 16 6.854 2 0.098 18 6.948 4 1.386 20 7.025 6 3.577 22 7.211 8 4.950 24 7.159 10 5.794 26 7.087 _ 12 . 6.317 28 7.203
14 6.671 30 7.242 As discussed in the text, the heat exchanger can be modelled as a ﬁrstorder system
with time constant 7'. Also, a time delay Td is introduced into the measurements
because the thermocouple sensor ‘is located some distance downstream of the heat
exchanger. Your job here is to ﬁt the model step—function response W) iET‘i) a[1——exp ] ,t>Td 0 , tSTd 1 to the sampled data using a sum—of—squared—errors objective function and the Nelder
Mead simplex optimizer discussed in lectures. The three parameters that are to be adjusted to the sampled data are (a, 7', Td). The ﬁtted response should look as shown
in the following ﬁgure. Heat Exchanger Step Response y, degrees Sampled Data In
Fitted Response 0 5 1O 15 20 25 30
time, minutes 4. A uniform cable, ﬁxed at two points and hanging under its own weight, takes up
the shape of a catenary. Starting at the lowest point, :1: = 0, the length of the cable for a position with positive :t can be expressed as % = sinhf and its vertical position as 3% = coshﬁ
where v is the length scale. The length of the cable is L = 100m and the horizontal distance between points
A and B is W 2 90m. Given that the the heights of points A and B are 40111
and 60 m respectively, determine the height of the lowest point of the cable. Also, determine the horizontal distance from point A to the lowest point of the cable. Hint: Set up the cable equations as a pair of nonlinear equations with unknowns v
and 13A. These equations will have the form f1('v, ma) = 0 and f2('v,1:a) = 0. Suitable values of the unknowns can then be determined by minimizing the objective function fobj = lf1l+lf2~ # hot_wire.py # Example of fitting hot—wire calibration data via linear least—squares.
# Exercise Sheet 8, Q1 # P.J. 12—May—2004 # from.matrix impdrt *
from math import sqrt print ’Begin hot4ﬂire exercise...’ U = [91.08, 85.04, 78.54, 68.51, 58.35, 49.63, 41.45, 34.32]
E = [0.8560, 0.8414, 0.8375, 0.8160, 0.7977, 0.7762, 0.7535, 0.7347]
na = 3 Ndata = len(U) # Set up the design matrix
phi = zeros(Ndata, na)
for i in range(Ndata): phi[i][0] = 1.0
phi[i][l] = U[i]
phiril[2] = sqrt(U[i]) # Set up normal matrix equation # Note that we need to build a list—of—lists from the original list E.
norm_mat = mmult(transpose(phi), phi) RHS = mmult(transpose(phi), transpose([E])) print 'normal matrix=’, norm‘mat print ’RHS=’, RHS # Solve the linear system to get the model parameters.
a = gauss_jordan(norm_mat, RHS)
print 'a=’, a ulist = [1; elist = [3
for i in range(35,9l):
u = float(i)
# Remember that a is a list—of—lists representing a column vector.
e = a[0][0] + aElJEOJ * u + a[2][0] * sqrt(u)
ulist.append(u); elist.append(e) import Gnuplot grph = Gnuplot.Gnuplot() d1 = Gnuplot.Data(U,E, with=’points’, title=’data’) d2 = Gnuplot.Data(ulist,elist, with=’lines’, title=’model’)
grph.title('Hot—wire calibration’) grph.xlabel('U, m/s') grph.ylabel(’E, volts') grph(’set xrange [30.0:90.0]’) grph(’set key bottom right’) grph.plot(dl, d2) grph.hardcopy("hot_wire.eps", terminal="postscript", mode="eps", fontsize=20) print ’Done.’ Begin hot_wire exercise... (with a little manual editing) normal matrix: [[8.0, 506.92000000000002, 62.884722431033026],
[506.92000000000002, 35153.343999999997, 4179.825340709589],
[62.884722431033026, 4179.825340709589, 506.92000000000002]] RES: [[6.4130000000000003], [412.71462600000001], [50.82160612049980811 a= [[0.47l99234251644046],
[—0.0012526560054564478], [0.052032632341551512]] Done. Hotwire calibration
0.86 0.84 0.82 0.8 E, volts 0.78 0.76 0.74 0.72 U, m/s #! /bin/env python # file: linesearch.py
'I'Ill Implements the bracket and golden_section line searches as described in Schiller and Harris Section 6.3. P.J.
08—Jan—04 from math import sqrt, pow def bracket(F, h=1.0, h_min=0.001, h_max=1000.0): Computes the Assumes F(x) This code is
In“ interval [a,c] that contains a minimum of f(x). decreases over range 0 <= x < x* and increases for x > x*.
a transcription of text of algorithm 6.3.1. .0
(F(b) >= F(a)) and (b > h_min):
i/.__. C while (F(c) < F(b)) and (c < hhmax): a = 0.0
b = h
c = 0
while
c b
b 2
.if c > 0.0:
return a,
else:
c = 2 * b
a — b
b = c c = 2 * c return a,c def golden_section_search(F, a, c, eps=l.Oe—8): Applies the Golden—Section search to find min(f(x)). This function is an almost direct transcription of algorithm.6.3.2.
errors in the original text were fixed.) (A couple of g = (3 _ sqrt(5))/2 bl = a + g*(c — a)
b2 = a + (l — g)*(c — a)
F1 = F(bl)
F2 = F(b2)
while (c — a) > eps*(a + c):
if Fl <= F2:
c = b2
b2 = b1
F2 = F1
bl = a + g*(c — a) # use c not b2
F1 = F(b1)
else:
a = bl
bl = b2
Fl = F2
b2 = a + (l — g)*(c — a) # use a not b2
F2 = F(b2) return (a + c)/2.0 # ________________ if __name__ == ’
print "Begin N = 0.6 main__':
test of linesearch.py..." def test_fun_l(x, N=N): diff = N — (x * x * x)
return diff * diff a,c= print
alpha
print
print print bracket(test_fun_l) “Bracket range: a=', a, ' c=', c
= golden_section_search(test_fun_l, a, c)
"After search, alpha=", alpha, ' F(alpha)=“, test_fun_l(alpha) “Expected alpha=", pow(N, “Done.‘ Begin test of linesearch.py...
Bracket range: a= 0.0 c= 2.0
After search, alpha= 0.843432662206 Expected alpha= 0.843432665302
Done.  1.0/3) F(alpha)= 4.36455681143e—17 #! lbin/env python
# file: schilling;ex691.py Solution to Exercise 6.9.1 in Schilling and Harris: Parameter identification for a heat exchanger. P.J. 08—Jan—04 II“! from math import exp
from.nelmin import minimize # sampled time stored as global variables
tk= [ 0.0, 2.0, 4.0, 6.0, 8.0, 10.0,
16.0, 18.0, 20.0, 22.0, 24.0, 26.0, 6.854, 6.948, 7.025, 7.211, 7.159, 7.087, def y(t,x):
"Model equation for heat exchanger response.“
a, tau, Td = x
if t > Td:
return a * (l — exp(—(t — Td)/tau))
else:
return 0.0 def least_squares(x): "Returns difference between sampled response and model response.“ global tk, yk sum = 0.0 for k in range(len(tk)):
diff = y(tk[k],x) — yk[k]
sum += diff * diff return sum if __name__ == ’__main__’:
print "Begin...“ 12.0, 14.0,
28.0, 30.0]
yk= [—0.036, e0.098, 1.386, 3.577, 4.950, 5.794, 6.317, 6.671,
7.203, 7.242] x, fx, conv_flag, nfe, nres = minimize(least_squares, [3.0,
[0.1, 3.0,
0.1, 3.0],
0.1]: 1.0e—9, 400) print “x=', x print “min—objefun—value=', fx print “convergence—flag=", conv_flag
print "number—of~fn—evaluations=', nfe
print “numberof—restarts=", nres
print ' —————————————————————————— ' # Generate a set of finely—spaced points to display the fitted response tp = [1; yp = [1
np = 10 * len(tk) dt = tk[—1]/np
for i in range(np):
t = dt * i tp.append(t); yp.append(y(t,x)) import Gnuplot
91 = Gnuplot.Gnuplot()
g1.title('Heat Exchanger Step Response") d1 = Gnuplot.Data(tk, yk, title=“Sampled Data",
d2 = Gnuplot.Data(tp, yp, title=“Fitted Response", with='lines 1") g1("set key bottom right")
g1.xlabel('time, minutes')
g1.ylabel(“y, degrees")
gl.plot(d1, d2) with="points 4") g1.hardcopy("schilling_ex691.eps", terminal="postscript", mode='eps‘, fontsize=20)
print l'Done." Begin... x= [7.2050186895088446, 4.2408048184998037, 3.0915626653490378]
minobj—fun—value= 0.0334851376832 convergence—flag= l number—of—fn—evaluations= 414 number—ofrestarts= 0 Done
Heat Exchanger Step Response
3
2
8’
'0
>2
Sampled Data a
. Fitted Response
0 5 10 15 20 25 30
time, minutes
~ t —T
mead, 936:) = CL ['1 ~ m4? ( d3 ] ‘t>Td
”C
2 o 0 \< Jc STd at Pam/rmi'm 6t, T Tcl I b . 0n~4
9 Tmmmlzmg €52) : Z [ 90% ) a?) N 9“]
k=® Begin catenary.py... z= [62.501793994386517, —32.328851670557292]
minobj~fun~value= 2.62900812231e12
convergence—flag= 1 number—of—fn—evaluations= 350
number—of—restarts= 0 Solution: Horizontal distance from A to low point: 32.3288516706
Mininmm height above ground: 31.4509210007
Height from reference to ground= 31.0508729937
Done. L .
.3: ~5mk EA } EL, scat»
49 0 1}
Ha
‘56 A = 005k 34% ; 99* “a = cask
U ’U 50!“ gr ’0, ac, ﬂu}: _._— #! /bin/env python
# file: caternary.py Solution to the catenary problem: Given (L, W, hA, hB) find the minimnm_height above the ground h_min
and the position of that low point with respect to point A (xA). P.J, 07~Jan—O4_ from.math import sinh, cosh
from nelmin import minimize L = 100.0 # specified parameters stored as global variables
hA = 40.0
hB = 60.0
W = 90.0 def objfun(z):
global L, hA, hB, W v = z[0] # rename independent vars to suit written formulation xA = z[1] # Now, set up the equations to solve. x1 = (W + xA) / v x2 xA / v _ fl v * (sinh(x1) — sinh(x2)) — L f2 = v * (cosh(xl) — cosh(x2)) — hB + hA # Combine the two equations as a penalty function.
f_obj = abs(fl) + abs(f2) return f_pbj I! if __name__ == ’__main__’:
print “Begin catenary.py...“
z, fz, conv_flag, nfe, nres = minimize(objfun, [60.0, 20.0], [1.0. 10].
1.0e—9, 400)
print "z=", 2
print 'min—obj—fun—value=', fz
print “convergence—flag=", conv_flag
print “number—of—fn—evaluations=', nfe
print 'number—of—restarts=", nres
print ' —————————————————————————— '
print “Solution:"
xA = z[1] print “Horizontal distance from.A.to low point:', —xA v = z[0] y_g = —hA + v * cosh(xA/v) h_min = v — y_g print "Minimum height above ground:“, h_min
print "Height from reference to groun =', y_g
print "Done." ...
View
Full Document
 Three '11
 PeterJacobs

Click to edit the document details