% Daniel Duhaney % Feb 2012 % IQP Tunnel Cost Modeling Octave script % This script is to be run in the same workspace that contains polyfitn.m and it's related functions % Tunnel data % face area f=[3.45967891,8.042477193,9.23,9.621127502,11.04466167,12.56637061,13.20254313,30.9156876,34.211944,34.211944,34.211944,40.71504079,52.8,60.13204689,68,85.62410777,86.29354105,108.8247695,109.3588403,128.6796351,157.5,208.25,216,264,284.75]'; % length l=[12500,2940,6475,3500,1720,2700,3000,6948,5600,5700,17000,25000,8500,10400,600,3400,9000,150000,9700,9600,7200,1220,334,2300,1066]'; % depth d=[39.5,70,27.4,20,35,17.5,30,21,30,30,30,20,53,17.5,10.75,65,30,45,30,60,35,15,2,24,22.6]'; % density (rho) r=[2067.5,2450,2050,2365,2265,2615,2450,1520,1806.7,1806.7,1806.7,2365,2625,1200,1520,1907.5,1675,1675,1550,2300,2365,1835,1200,2615,1333.3]'; % outturn cost c=[527401421.3,106219588.5,285000000,364070520.3,58827057.12,160386268.6,113100795.5,1877777778,1260316074,978556.6729,2043451963,2451439030,856550238.2,1233442623,530000000,3375646567,860812372.2,13251396648,727671294.8,2925333753,950285544,186650238.5,569698449.6,909081273.9,232630894.3]'; % cost/m^3 cm3=[12195.38426,4492.286503,4768.737163,10811.6381,3096.679501,4727.086566,2855.530544,8741.892776,6578.300307,5.018031401,3513.48221,2408.386663,1908.534399,1972.330202,12990.19608,11595.3001,1108.377498,811.7880213,685.9770922,2368.069087,837.9943068,734.6554561,7896.685097,1497.169423,766.3840414]'; % inputs for test tunnels testf=[6.157521601,30.9156876,32.16990877,52.16810951,95.03317777,108.8247695,358.1]; testl=[3300,6948,1249,15000,5600,150000,680]; testd=[21,22.1,15,21,34,45,9]; testr=[1600,1573.3,1935,1760,1985,1825,1200]; testcm3=[7486.528961,8741.892776,29865.44419,546.5202553,481.8556235,811.7880213,3520.946099]; % inputs for 'prediction' tunnel set predf=[108.618681,9.621127502,35.36184545]; predl=[5899,9000,7242]; predd=[192.5,80,36.6]; predr=[2800,2700,2650]; % evaluates and stores and presents data on polynomials up to degree 12 for i=1:8 disp('DEGREE'),disp(i) model5D(i)=polyfitn([f,l,d,r],c,i); %approximates 5Dimension dataset with polynomial of degree i disp('5D POLYNOMIAL') disp('RMSE:'), disp(model5D(i).RMSE) % displays RMSE for 5d polynomial of degree i disp('') model3D(i)=polyfitn([d,r],cm3,i); %approximates 3dimensional dataset with polynomial of degree i disp('3D POLYNOMIAL') disp('RMSE:'), disp(model3D(i).RMSE) % displays RMSE for 5d polynomial of degree i disp('') disp('_____________') disp('') endfor % to access the coefficients and model terms of an individual model of degree i, enter model5D(i).Coefficients and model5D(i).ModelTerms, respectively. % see polyfitn.m for more details % evaluates model output for tunnels in the dataset for j=1:25 m(j)=polyvaln(model5D(3),[f(j) l(j) d(j) r(j)]); n(j)=polyvaln(model3D(5),[d(j) r(j)]); o(j)=polyvaln(model3D(3),[d(j) r(j)]); endfor % evaluates model output for test tunnels not in the dataset for k=1:7 s(k)=polyvaln(model5D(3),[testf(k) testl(k) testd(k) testr(k)]); t(k)=polyvaln(model3D(5),[testd(k) testr(k)]); u(k)=polyvaln(model3D(3),[testd(k) testr(k)]); endfor % evaluates model output for data outturn cost imputation for w=1:3 x(w)=polyvaln(model5D(3),[predf(w) predl(w) predd(w) predr(w)]); y(w)=polyvaln(model3D(5),[predd(w) predr(w)]); z(w)=polyvaln(model3D(3),[predd(w) predr(w)]); endfor % function handles representing polynomials % C3(d,r) c3=@(d,r) (1.60785249816796e-001)*(d.^3)+(-9.67853886666806e-003)*(d.^2).*r+(6.18204392443470e+000)*(d.^2)+(-2.98502603295617e-004)*d.*(r.^2)+(1.92029273329062e+000)*d.*r+(-2.39236650599325e+003)*d+(8.87468360175018e-006)*(r.^3)+(-5.28650499314227e-002)*(r.^2)+(9.16445204287247e+001)*r+(-3.81802502457303e+004); % C5(d,r) c5=@(d,r) (7.12033097624903e-003)*(d.^5)+(1.99441192088892e-004)*(d.^4).*r+(-1.83101314495651e+000)*(d.^4)+(-3.91919219461148e-005)*(d.^3).*(r.^2)+(1.41487331356153e-001)*(d.^3).*r+(-1.62444332151579e+001)*(d.^3)+(8.28224229444261e-007)*(d.^2).*(r.^3)+(-1.71096427039142e-003)*(d.^2).*(r.^2)+(-2.96576494690989e+000)*(d.^2).*r+(2.17445571725495e+003)*(d.^2)+(5.28560972369608e-009)*d.*(r.^4)+(-7.71000033594450e-005)*d.*(r.^3)+(2.44901711861428e-001)*d.*(r.^2)+(-2.14442923306664e+002)*d.*r+(4.90406404484917e+004)*d+(1.68493005361992e-010)*(r.^5)+(-1.80361479258466e-006)*(r.^4)+(7.95473848195651e-003)*(r.^3)+(-1.70936117847365e+001)*(r.^2)+(1.70178164672124e+004)*r+(-6.25905106749705e+006); % O(f,l,d,r) O=@(f,l,d,r) (-2.96255941113289e+004)*(f.^3)+(-4.44377174708671e+001)*(f.^2).*l+(-9.06755894474798e+005)*(f.^2).*d+(3.44425203091532e+004)*(f.^2).*r+(-2.39173190952879e+007)*(f.^2)+(-3.08730662110037e-001)*f.*(l.^2)+(6.42937166532731e+002)*f.*l.*d+(-7.85808917955927e+001)*f.*l.*r+(-1.52320113557384e+004)*f.*l+(-4.60783918046207e+006)*f.*(d.^2)+(2.75032860951945e+005)*f.*d.*r+(7.66481875122618e+006)*f.*d+(-3.50626095764708e+003)*f.*(r.^2)+(-1.84143646625962e+006)*f.*r+(4.29347531122274e+009)*f+(-7.35181676227068e-003)*(l.^3)+(1.37228425682328e+001)*(l.^2).*d+(4.87248418192956e-001)*(l.^2).*r+(4.71589798292492e+001)*(l.^2)+(-3.69517689718922e+004)*l.*(d.^2)+(1.24232966009547e+003)*l.*d.*r+(-3.17929064791201e+005)*l.*d+(2.04339048921672e+000)*l.*(r.^2)+(-2.49370453702774e+004)*l.*r+(-8.28926328654429e+006)*l+(7.75855830346090e+006)*(d.^3)+(-1.35843214078337e+004)*(d.^2).*r+(-7.55825606798212e+008)*(d.^2)+(-9.21237973834454e+003)*d.*(r.^2)+(3.44696102332892e+007)*d.*r+(-1.11607836988330e+010)*d+(1.46407991645023e+002)*(r.^3)+(-5.28402897761373e+005)*(r.^2)+(3.93425458953762e+008)*r+(9.43062442108692e+010); % create meshgrid for mesh plot % see close-up of mesh plots, limit the length of the linspace lin1 to linspace(10,30,200)' lin1=linspace(10,35,200)'; lin2=linspace(1190,2635,200)'; [depth,density]=meshgrid(lin1,lin2); % calculate cost with one of the C polynomials cost=c3(depth,density); % plot mesh mesh(depth,density,cost); hold on % includes scatter plot points for dataset tunnels plot3(d,r,cm3,'g.'); hold on % includes scatter plot points for test tunnels % plot3(testd,testr,testcm3,'r.'); % labels grid xlabel('depth (m)') ylabel('substrate density (kg/m^3)') zlabel('cost per metre-cubed (US$)') title('Scatter Plot and Data fit of Empirical Tunnel Cost Data')