diff --git a/README.md b/README.md index ed66e19..8437f08 100644 --- a/README.md +++ b/README.md @@ -86,4 +86,27 @@ Planned parameters include: - flange width - local reinforcement - mesh density -- joint and spot-weld idealization \ No newline at end of file +- joint and spot-weld idealization + +## Stochastic Engineering + +BodySimPy includes Monte Carlo uncertainty propagation for selected manufacturing, material and loading parameters. + +The baseline stochastic study evaluates 500 independently sampled configurations with uncertainty in: + +- wall thickness +- Young's modulus +- applied load + +Each sampled configuration is evaluated using automated CalculiX static and modal analyses. + +Reported quantities include: + +- response mean and standard deviation +- 5th and 95th percentiles +- stress-threshold exceedance fraction +- mode-1 frequency variation +- correlation analysis +- standardized linear sensitivity coefficients + +The stochastic model is intended as an engineering uncertainty study for the simplified structural surrogate and does not represent production vehicle reliability. \ No newline at end of file diff --git a/configs/stochastic_crossmember.yaml b/configs/stochastic_crossmember.yaml new file mode 100644 index 0000000..cdd2471 --- /dev/null +++ b/configs/stochastic_crossmember.yaml @@ -0,0 +1,42 @@ +project: + name: stochastic_crossmember + +geometry: + type: rectangular_hollow_section + length_m: 1.0 + width_m: 0.080 + height_m: 0.040 + thickness_m: 0.0015 + +material: + youngs_modulus_pa: 2.10e11 + poisson_ratio: 0.30 + density_kg_m3: 7850 + +loading: + tip_force_n: 1000 + +mesh: + elements: 40 + +analysis: + static: true + modal: + modes: 10 + +stochastic: + samples: 500 + seed: 42 + stress_threshold_pa: 3.50e8 + + thickness_m: + mean: 0.0015 + standard_deviation: 0.00005 + + youngs_modulus_pa: + mean: 2.10e11 + standard_deviation: 5.0e9 + + tip_force_n: + mean: 1000.0 + standard_deviation: 80.0 \ No newline at end of file diff --git a/docs/figures/frequency_histogram.png b/docs/figures/frequency_histogram.png new file mode 100644 index 0000000..171872e Binary files /dev/null and b/docs/figures/frequency_histogram.png differ diff --git a/docs/figures/frequency_sensitivity_ranking.png b/docs/figures/frequency_sensitivity_ranking.png new file mode 100644 index 0000000..c4f47c1 Binary files /dev/null and b/docs/figures/frequency_sensitivity_ranking.png differ diff --git a/docs/figures/stochastic_correlation_matrix.png b/docs/figures/stochastic_correlation_matrix.png new file mode 100644 index 0000000..2a63360 Binary files /dev/null and b/docs/figures/stochastic_correlation_matrix.png differ diff --git a/docs/figures/stress_histogram.png b/docs/figures/stress_histogram.png new file mode 100644 index 0000000..e2389a9 Binary files /dev/null and b/docs/figures/stress_histogram.png differ diff --git a/docs/figures/stress_sensitivity_ranking.png b/docs/figures/stress_sensitivity_ranking.png new file mode 100644 index 0000000..91159fd Binary files /dev/null and b/docs/figures/stress_sensitivity_ranking.png differ diff --git a/docs/validation/stochastic_correlation.csv b/docs/validation/stochastic_correlation.csv new file mode 100644 index 0000000..8232404 --- /dev/null +++ b/docs/validation/stochastic_correlation.csv @@ -0,0 +1,6 @@ +,thickness_mm,youngs_modulus_gpa,tip_force_n,max_stress_mpa,mode_1_frequency_hz +thickness_mm,1.0,0.08519099399636426,-0.0006070039120654317,-0.3777326617507689,-0.025829604838853208 +youngs_modulus_gpa,0.08519099399636426,1.0,0.02851868542650945,-0.00796082089382417,0.993791474405028 +tip_force_n,-0.0006070039120654317,0.02851868542650945,1.0,0.9254327505427687,0.029427867937940846 +max_stress_mpa,-0.3777326617507689,-0.00796082089382417,0.9254327505427687,1.0,0.03472301484721318 +mode_1_frequency_hz,-0.025829604838853208,0.993791474405028,0.029427867937940846,0.03472301484721318,1.0 diff --git a/docs/validation/stochastic_samples.csv b/docs/validation/stochastic_samples.csv new file mode 100644 index 0000000..e02d8cf --- /dev/null +++ b/docs/validation/stochastic_samples.csv @@ -0,0 +1,501 @@ +sample,thickness_mm,youngs_modulus_gpa,tip_force_n,max_stress_mpa,mode_1_frequency_hz +1,1.5152358539877218,204.80007946879752,1060.0360956645166,198.608,48.66397 +2,1.5470282358195606,200.2448240567308,895.8256394510146,164.74,48.08207 +3,1.5063920201583645,208.41878703828206,998.6559073996569,188.095,49.1027 +4,1.457347803621321,214.39698987431416,1062.2233548343158,206.128,49.86205 +5,1.5033015348780607,215.63620603484014,1037.4007473801637,195.754,49.94947 +6,1.4570353768558382,211.8437539204125,923.2893919336801,179.202,49.56464 +7,1.5439225150653637,209.75037044506874,985.2110109163791,181.505,49.21382 +8,1.465953522779803,216.11270669337014,987.6376414344958,190.638,50.05058 +9,1.4785836088918447,208.23933224755885,1042.5847348442678,199.693,49.11515 +10,1.5182722032182039,212.06366305797994,1034.4656802406307,193.468,49.51572 +11,1.607082380043523,207.96792491807693,959.020581674277,170.453,48.92802 +12,1.4593113635876063,213.07989711287746,1090.3177834176713,211.322,49.70626 +13,1.4943026271172561,205.79921761518736,934.0415027447009,177.206,48.80766 +14,1.532529639391235,213.71627085601722,1043.4523414644157,193.517,49.69084 +15,1.4667245146355652,211.1608066153336,1009.3348647312582,194.734,49.4729 +16,1.5109344298364507,214.35714388974094,1017.8876439019746,191.199,49.79175 +17,1.5339456781535947,210.33789534744446,1023.1295518951987,189.591,49.29481 +18,1.5315644112919269,202.71422090072167,974.4263026914159,180.818,48.39606 +19,1.4764813672853603,206.8056107587833,977.9886199018653,187.561,48.94831 +20,1.5747470655617197,205.6708444215338,1077.4622683673185,195.015,48.69588 +21,1.4158565114192099,208.32557485007115,1013.0202452084004,201.783,49.20109 +22,1.529311116567964,213.55613289896425,1063.467778815994,197.602,49.67616 +23,1.4825637463875783,207.68824103667714,1068.6380705005722,204.187,49.04531 +24,1.4904347837559193,203.6215683833104,909.3370228797215,172.922,48.55337 +25,1.4540273856999193,212.48580372026882,1011.3940588856452,196.669,49.64336 +26,1.5345242677033886,207.86373676831732,1012.6831752861372,187.592,49.00333 +27,1.5312795196983668,208.4532673013988,1036.542019004593,192.377,49.07669 +28,1.4669037029466676,208.18473076717464,969.4609684801337,187.02,49.12281 +29,1.440208017720548,212.43486240392792,962.4478127837821,188.773,49.65427 +30,1.5006247059363844,212.40373329452953,1035.7224940823955,195.751,49.57694 +31,1.5332692554486393,209.5075725774529,966.1361350364677,179.101,49.19823 +32,1.49601408945468,201.56332783020983,884.231002206153,167.583,48.30072 +33,1.43386501938228,205.0137658619926,1031.981938137875,203.222,48.78686 +34,1.454726047231997,208.10918722980304,1103.9382638228853,214.571,49.1286 +35,1.4821868014469286,213.68757784233543,925.3105855992098,176.842,49.7491 +36,1.4897281221066185,205.24988972544708,972.877353927955,185.084,48.74797 +37,1.54201540687287,201.36339788403828,1034.7538914836687,190.844,48.22213 +38,1.511886780116614,207.029250221516,884.3153716489236,166.015,48.93213 +39,1.5036064753856935,207.352536454681,1018.6140969083763,192.174,48.9803 +40,1.501092607276172,218.0088944566046,980.851549802158,185.329,50.22625 +41,1.4488251253689068,210.89637817478157,1017.5997347177412,198.517,49.46367 +42,1.5679593787620218,214.17555622957292,1028.5496847319607,186.884,49.70081 +43,1.573165144560978,204.05618472838574,948.8198773800202,171.886,48.50624 +44,1.4536712029297236,208.05095098422117,889.8651081954953,173.075,49.123 +45,1.5317575473407201,208.88888651450614,882.3354964397873,163.711,49.12737 +46,1.4492210459396229,211.56756923725098,1067.0501254315504,208.113,49.54183 +47,1.5998365445845892,224.56931233003664,1033.1527546620798,184.371,50.8525 +48,1.4505230939984066,199.33976859634345,1021.4169169875069,199.051,48.08732 +49,1.4593529452344838,207.92321369910155,951.0322560752154,184.322,49.10106 +50,1.492960455679181,215.32990115393824,1012.5638853956276,192.259,49.92668 +51,1.4920682581480658,204.82173123587094,866.0253644236515,164.524,48.69429 +52,1.4756846045463334,209.73108724590838,1141.4343930863906,219.013,49.29427 +53,1.5065137260736443,214.9136975551153,960.0563521168679,180.812,49.86177 +54,1.4407528116791488,205.17441618883814,941.9819148371397,184.696,48.79771 +55,1.6064234866217584,205.89306660387808,1067.0791362989075,189.728,48.68413 +56,1.4548536410956487,214.6578650643712,1030.796077288469,200.339,49.89544 +57,1.4921681051170954,209.7961873693228,947.6169843656488,180.013,49.28194 +58,1.522303610074104,207.72508259829613,901.9515388986201,168.284,49.00172 +59,1.436103121284019,210.8629395886106,1126.3273005128349,221.488,49.47521 +60,1.507999580678672,209.4068083694506,1022.8660911682034,192.47,49.21701 +61,1.5653000870853413,211.09691250681928,967.1258215333011,175.99,49.34554 +62,1.5553144355029944,212.14378219230807,1122.8604793596794,205.505,49.47992 +63,1.5091617218610953,203.8776548413975,890.5472640603468,167.456,48.56151 +64,1.5825463966115625,218.61832860391647,985.638462937392,177.609,50.19559 +65,1.4808406339432008,217.307221461211,911.4363454365209,174.333,50.17034 +66,1.4552636490522086,213.21663397344523,968.4315901712329,188.17,49.72715 +67,1.4997439066639964,209.18278550737745,1027.0059639039148,194.207,49.20067 +68,1.5703740930656858,210.45292453450872,1051.5151034621485,190.793,49.26405 +69,1.397491394948449,209.7564079903494,932.5415783765703,187.965,49.39202 +70,1.4390593469788187,205.60923816535623,973.2701247439351,191.034,48.85143 +71,1.5457951271178005,203.36803641130217,1002.4505194075534,184.48,48.45705 +72,1.475791528333321,208.3616345281902,1080.2206260243684,207.254,49.13294 +73,1.5269057718501962,216.68699053721372,987.6395456600079,183.772,50.04194 +74,1.465202869416465,208.8807059155975,1019.3997433016978,196.86,49.20691 +75,1.5088286679226854,204.57805963883317,1007.23918253023,189.436,48.64525 +76,1.5114114165069452,222.58737018766962,1150.1475689025335,215.981,50.73803 +77,1.4573378324720592,208.5630831922541,882.9246398530397,171.335,49.17899 +78,1.4704646493018256,211.5780250179517,1096.4682896710588,211.059,49.5172 +79,1.4635458081128196,206.72926779966102,828.2168776209076,160.103,48.95484 +80,1.4918667039727544,204.68792794070217,957.6448458107141,181.951,48.67862 +81,1.4561569610916207,209.52868722873723,859.3817286914696,166.889,49.29413 +82,1.4266477377304505,220.6462355601415,897.0061934980776,177.45,50.62169 +83,1.4451607210772681,219.18456764160658,1232.4053735392326,240.973,50.4308 +84,1.4414216685587329,208.15875521605986,1027.3244440752408,201.345,49.15052 +85,1.5864348822202798,205.06571460785884,980.3777723246312,176.273,48.61018 +86,1.5388668788030873,212.1738303723331,969.9075143015921,179.213,49.50347 +87,1.493308851774198,203.12552095815008,980.9461004820267,186.216,48.49076 +88,1.4866806254995524,211.16084944813127,955.5738224944787,182.128,49.44863 +89,1.523576926127257,215.06357908909916,1012.4343462134773,188.755,49.85821 +90,1.5175878204199602,210.26577673788935,1000.0067514473131,187.1,49.30621 +91,1.4639220983228576,211.58247130836654,992.2170721269208,191.762,49.52568 +92,1.6046584154465298,217.86677451237622,1030.86772420452,183.469,50.08192 +93,1.4618471395152617,204.43794264008292,1095.2914362471092,211.954,48.68481 +94,1.513137461257357,212.40071701695805,860.4331210420646,161.411,49.56132 +95,1.546371924077909,212.2721016910718,911.1655452468417,167.625,49.50578 +96,1.4764237596275025,211.31858602178258,1004.1973438684995,192.594,49.47958 +97,1.4853914407098223,209.4825586596702,979.8418097434492,186.9,49.2533 +98,1.5076281256051236,217.3574598649658,794.6673247254962,149.564,50.14308 +99,1.4881574867745158,210.88256210686225,1023.6795191749679,194.934,49.41424 +100,1.4814042709339357,201.2163910876071,1026.2396386971288,196.225,48.27648 +101,1.5863675107082091,202.3306929754193,1069.106241095855,192.234,48.285 +102,1.4835737388438555,209.69337827085576,915.7681191436841,174.871,49.28028 +103,1.4832771913822063,216.50022297349477,1046.6124193590165,199.892,50.0741 +104,1.5866155802704973,215.88705944488845,1035.126934314366,186.098,49.87605 +105,1.5871967263083897,212.1949657936915,1066.2390541437646,191.629,49.44701 +106,1.4851714523232646,210.33272907929603,944.2060941359779,180.126,49.35341 +107,1.5494791967297425,204.10848188442566,1062.5880273934129,195.131,48.54078 +108,1.490467447124742,215.8562354673366,1060.0695190985007,201.582,49.99072 +109,1.59103230788234,213.65387345569323,874.2367795477408,156.783,49.61201 +110,1.49665234135509,204.13996403702592,958.537612610698,181.596,48.6077 +111,1.575561421663836,213.18766905448848,944.0855657137002,170.796,49.57676 +112,1.44931413167798,210.16391046438133,902.6751880495098,176.044,49.3771 +113,1.4664429861301709,211.56004735775622,1092.4249628885357,210.801,49.51999 +114,1.530438073611222,198.5435524570095,1024.3493383224209,190.207,47.89694 +115,1.5036016786616289,212.06945142686345,1129.2967745320966,213.056,49.53428 +116,1.396838086159763,207.04448287404446,1047.2725058234685,211.179,49.07248 +117,1.420920280027303,217.37974524026754,1029.4685290160526,204.398,50.25265 +118,1.5423291992921704,207.14528169239102,1065.101095197297,196.405,48.90917 +119,1.5534235777994518,211.16439010068586,1018.7520714158957,186.655,49.36787 +120,1.5135171619640864,205.6832736763031,988.1977055807718,185.337,48.77085 +121,1.4923738733773821,211.91696932176748,1079.9859397528235,205.133,49.53015 +122,1.4470731959010796,209.3749548484021,1118.5164438127179,218.443,49.28705 +123,1.462820588559372,205.8887499135484,1016.1844953464955,196.527,48.85609 +124,1.542219259524266,210.05713027780328,1106.3168472349591,204.018,49.25186 +125,1.5428396991283264,214.2091003344503,1044.3293201059275,192.518,49.73547 +126,1.6163826550173042,208.97419155566976,839.7182166346789,148.482,49.035 +127,1.5802127180758572,207.71150284521588,1008.6304355474426,181.992,48.93027 +128,1.5654775342430263,201.98870229470538,899.8682228715088,163.735,48.26903 +129,1.4199361031440574,206.0293185498729,1035.1709284813594,205.659,48.92428 +130,1.5262093921406386,211.38137089664275,886.9787292886107,165.11,49.42634 +131,1.3844948281683678,210.2717679269519,962.2579173089331,195.606,49.46843 +132,1.5229692885858306,213.50976813110634,1011.0593145646945,188.567,49.67852 +133,1.5380066542687927,211.14605687057534,1042.4051764481173,192.705,49.38448 +134,1.464766336856506,209.1019429308427,1015.7420877325566,196.206,49.23349 +135,1.541026423770871,208.03129413887714,1041.693380458569,192.234,49.01523 +136,1.486708060404231,209.41228916335933,1066.3615233619257,203.24,49.24344 +137,1.4003469814472922,203.5176383596262,881.4251682045743,177.333,48.64854 +138,1.3833191940075849,206.60867779922413,1059.9547119782192,215.632,49.03706 +139,1.4857557966808712,210.98895040970925,1087.1373997368635,207.32,49.42963 +140,1.5663843066133845,209.65431032636928,1108.2868711655517,201.553,49.17532 +141,1.5046063329217054,205.8130088808627,952.4479718241012,179.584,48.79693 +142,1.4259731743717492,205.55933073133173,971.3586649540457,192.241,48.8612 +143,1.54017925096893,218.60384915582992,889.4254787698368,164.216,50.24636 +144,1.519641373412996,204.79728030421234,1037.9757670770562,193.968,48.65836 +145,1.4934456674656138,200.84547087076237,1074.2637593254788,203.914,48.21768 +146,1.4697499643595635,207.33049880813223,914.4198070968248,176.094,49.0185 +147,1.4672858361656222,212.13945220347458,984.8604527250876,189.946,49.58672 +148,1.516433100114124,211.80960926964423,1105.6529324422534,207.007,49.48829 +149,1.482860692456781,202.61571091577133,1085.3777933257586,207.348,48.44232 +150,1.4834255913951373,215.57296222288687,1030.6701690545976,196.83,49.96657 +151,1.4934431234898329,211.74387947023146,1156.0810048101039,219.445,49.50862 +152,1.6038490264526877,210.3469056756639,1012.815247465366,180.338,49.211 +153,1.5538120078733193,205.77169483616328,1026.6456298094604,188.059,48.73295 +154,1.498706857602307,211.56954105728778,933.3304955240476,176.602,49.48182 +155,1.4205216253345516,199.63508282004355,910.6092709608314,180.844,48.15841 +156,1.4770662357530333,208.5340420667562,1154.9784929943614,221.425,49.15172 +157,1.555299668495365,205.18954441829186,1027.8166761960765,188.112,48.66218 +158,1.4796460881824838,208.57818080979953,1014.8260519532306,194.25,49.15381 +159,1.5309585558487695,208.30370758057992,1085.1081226187487,201.428,49.05947 +160,1.4429030886928798,210.03169531181135,1207.8138981276666,236.499,49.36934 +161,1.5111539871312638,217.16607253303084,1007.3216142538566,189.19,50.11665 +162,1.52903885476649,209.71608402803653,986.3673935068637,183.305,49.22784 +163,1.461025880163729,212.1515067947878,931.8770248651385,180.422,49.59576 +164,1.5332792618156714,215.42643502222955,1029.3225125789784,190.813,49.88834 +165,1.4856875633222153,212.2698278965435,975.306155556287,186.001,49.57952 +166,1.5467773562732574,200.84296957892437,973.1514105505083,178.987,48.15412 +167,1.4004594002438,202.52469584886398,1109.1089783851128,223.124,48.52958 +168,1.5447592490985844,206.40259883357604,879.7997235167282,162.007,48.81849 +169,1.3517735581079173,207.28252246033682,1193.6332009797923,247.974,49.15497 +170,1.5217442135731825,207.20213856975258,1037.2064167600245,193.583,48.94068 +171,1.4219520823502778,208.51338318651182,1007.9581978412591,199.995,49.21589 +172,1.4956949674089475,213.9540306084504,1027.5716180988843,194.787,49.76357 +173,1.5334163009053998,206.5581385888462,1071.8252326728439,198.676,48.85052 +174,1.5814468473811996,205.14925240174296,928.9843475428352,167.504,48.62606 +175,1.5667892168160116,209.0432800665247,1112.3057113653244,202.237,49.10312 +176,1.4778732144053908,217.27522788135354,1010.5188653444362,193.635,50.17031 +177,1.5129114411661346,217.82359010834952,971.0583618044664,182.188,50.1903 +178,1.4529438952037521,207.7571789598582,1036.1867160514512,201.626,49.08918 +179,1.4217120463909743,213.18735451325563,956.8982945905806,189.893,49.76474 +180,1.5573906330366782,198.02869847549766,937.0747379916666,171.298,47.80306 +181,1.4156765924382948,205.86885266818024,1019.8132720887154,203.16,48.91034 +182,1.4910386687275123,208.73311215525993,987.2652102895952,187.673,49.15828 +183,1.5101694120309972,204.95731979028446,1056.547971271924,198.551,48.68872 +184,1.5331332985192743,211.9251896882805,1044.522675420101,193.647,49.48145 +185,1.5148209000404331,220.17536651365535,993.0324663157992,186.1,50.45814 +186,1.484645839083271,206.23236209821312,917.4189863537731,175.072,48.87061 +187,1.4377764106199462,205.5560134338454,994.3455694678343,195.328,48.84664 +188,1.5167147564248857,210.2557102927608,938.757177805623,175.731,49.30609 +189,1.5450092282009982,213.69706336150472,987.2281354062386,181.762,49.67334 +190,1.4673541927667646,212.74213960414897,1015.0378869988883,195.758,49.65703 +191,1.4275936370292477,209.66009872007797,1020.9628649659505,201.851,49.3442 +192,1.4550152606756124,210.94921696418723,883.6142011793769,171.716,49.46234 +193,1.5668093050006084,216.23974992529716,979.7986132557628,178.143,49.94115 +194,1.5181727168919537,197.9503901710006,907.4921871787734,169.731,47.83979 +195,1.4853110539923935,204.6393348928657,1057.1517186104527,201.656,48.6807 +196,1.5998648265373998,204.1169264028535,933.0029276731846,166.496,48.48153 +197,1.5117724184154964,218.0555807424981,902.2100549968077,169.386,50.21842 +198,1.5124518061153471,219.10649425406555,867.8592681456587,162.871,50.33845 +199,1.4359465396577085,207.8819669676958,958.3529269715671,188.475,49.12443 +200,1.5406300643876822,211.20829859910418,858.0030352286298,158.373,49.38856 +201,1.5257705201450407,207.11230556322175,1101.9557773552992,205.18,48.92522 +202,1.468620623178154,206.81692358448882,1043.2905288643626,201.05,48.95911 +203,1.538146324116847,212.24049681968103,865.1522146092153,159.924,49.51213 +204,1.5269017219964125,204.828459741628,1018.822088927061,189.575,48.65337 +205,1.428813278230047,212.23161073771533,935.4720876554102,184.806,49.6444 +206,1.4358682674869738,213.569100682412,1019.3315617341401,200.477,49.79196 +207,1.4693011599303551,217.2558942451053,964.7478062744206,185.837,50.17865 +208,1.5016053836985468,211.344567238393,950.4267252697109,179.525,49.45198 +209,1.5235568146699965,207.33273826417638,967.0689342227033,180.3,48.95392 +210,1.568132131990167,204.79706973743782,806.9775735256408,146.61,48.60025 +211,1.5805468497505948,222.746639763035,967.5784587640535,174.552,50.66981 +212,1.4031580979689908,208.4475801190203,977.1021640127465,196.225,49.23084 +213,1.4905038079488873,204.43305979039076,1046.3648913981183,198.972,48.64994 +214,1.5262253687643044,202.52797190090348,1055.935738535591,196.559,48.38018 +215,1.6026342490972947,210.8598016621853,973.0139870345195,173.368,49.27244 +216,1.492899839279624,213.0762838348428,861.5462715594225,163.591,49.6648 +217,1.5082195351488652,208.04768025268595,1147.8260010364836,215.956,49.05676 +218,1.4912913633097522,218.33943805559758,911.7007441774394,173.282,50.27643 +219,1.5293629583401398,211.59700131895616,930.4762201708174,172.886,49.44771 +220,1.508869805731847,216.062594184158,974.096663657072,183.199,49.99197 +221,1.4154018675415119,209.91218586761275,927.8061524165989,184.863,49.38864 +222,1.482882953066806,209.5920611515744,863.5478227187587,164.968,49.26921 +223,1.419217083033966,212.41033417762762,958.1825043064256,190.451,49.677 +224,1.3717627829089718,213.9242201830407,1021.7895804219391,209.458,49.9118 +225,1.4643062587981759,203.4158489447088,1066.864631705147,206.14,48.56003 +226,1.5174675311280215,221.91301141342038,1033.615087897995,193.401,50.65356 +227,1.519385157063101,209.1653603494572,1065.3420693760634,199.112,49.17484 +228,1.5312542600624077,216.258625042826,958.2941665090173,177.857,49.9871 +229,1.4782296263528019,207.6044841454401,1063.264137694319,203.698,49.04065 +230,1.5749187220148366,207.70579753425352,966.0181017765527,174.827,48.93598 +231,1.5157038611934899,208.7711924999355,1076.1642922511912,201.574,49.13294 +232,1.3874111354690477,205.86647476541552,937.4066932830843,190.191,48.94399 +233,1.3839821993806927,205.18180777815599,926.7875107323259,188.459,48.86665 +234,1.4899447672740584,215.56482945830132,980.3951479047438,186.49,49.95761 +235,1.4484593139864423,209.71522729393234,1083.933889715308,211.506,49.3254 +236,1.4512020591942547,205.44712558039276,1044.683916741938,203.499,48.8176 +237,1.48892319421346,213.23742345354736,998.9082566563762,190.129,49.68844 +238,1.5350831824044784,204.8246094297968,999.0332279369386,185.003,48.64311 +239,1.4894654992818421,203.92054537947934,874.921764170505,166.475,48.59016 +240,1.534287451763733,208.2452380889101,918.2175774682137,170.116,49.04856 +241,1.4951910531538335,215.64009428674652,817.5409723702462,155.021,49.95988 +242,1.4251680655048335,205.38556779052365,1116.8943093376263,221.157,48.84151 +243,1.5141293491459649,213.8365862023061,908.7871264699015,170.382,49.72735 +244,1.4440232062364737,212.23906858361354,1004.6619464748358,196.582,49.62673 +245,1.5274369409432065,209.06164504535772,1022.2514981163702,190.153,49.15291 +246,1.5079059542452555,213.88883694542463,1064.5606625555577,200.327,49.74104 +247,1.4190064002804115,198.76365722329572,1080.1396318121206,214.72,48.05497 +248,1.5593862546630968,204.89688490243333,851.2131674763498,155.424,48.62257 +249,1.50495174105036,214.65419113175267,1143.8075617947452,215.62,49.83358 +250,1.5258148786050725,208.14141651707152,928.549550273211,172.888,49.04657 +251,1.5005725645835328,208.50367648057372,918.7945675416947,173.657,49.11974 +252,1.6024377825431644,218.92584201036067,1090.8838946844344,194.391,50.20625 +253,1.4539574801095436,214.2750966214864,1051.1701127598321,204.412,49.85203 +254,1.5221272812102027,216.24832633521888,1050.8296914289022,196.081,49.99714 +255,1.5370007118152167,213.18453141686106,1027.2633126401686,190.017,49.62352 +256,1.410819443514291,210.4181053870528,955.50464329207,190.942,49.45369 +257,1.436007952117999,218.40908309650413,1138.3196272708562,223.859,50.35283 +258,1.5679610300147806,211.27606708427936,1108.0500096762125,201.328,49.36324 +259,1.5006026590300914,211.0139863844889,912.5222903630414,172.469,49.41451 +260,1.5198495652249422,210.30192964212208,895.7878306532428,167.376,49.3077 +261,1.4974401444091703,209.60135218071593,1143.8048933315713,216.593,49.25266 +262,1.5447106664256705,210.05722719545733,1019.9029847840866,187.81,49.24885 +263,1.50221061890521,208.98543007602225,913.405824208937,172.469,49.17447 +264,1.4924474061196535,206.26950871673742,899.9747568861483,170.934,48.86563 +265,1.5255610907273767,211.95632329966617,857.0633990905537,159.602,49.49431 +266,1.4938657681252436,214.97850692949962,1084.7378991117057,205.85,49.88481 +267,1.5512918421442619,210.1945660345073,932.3963737624679,171.043,49.25696 +268,1.445815058632097,211.7230245237518,1030.342402749213,201.381,49.56418 +269,1.5643665555802533,215.49989825171042,989.4217074775063,180.144,49.85864 +270,1.4377897060569564,208.40446701494548,1017.378495584842,199.851,49.1839 +271,1.4898958377305298,207.1105294596117,1020.2307851975919,194.073,48.96822 +272,1.4748021129294695,206.85966958137314,1024.9162267442753,196.762,48.95673 +273,1.4799003709935399,211.22052194489174,1021.8563218281698,195.565,49.46387 +274,1.4430285499454214,207.59379438067293,1115.022155193317,218.313,49.08183 +275,1.4418960243632857,199.41672775514735,851.0523877405991,166.748,48.1068 +276,1.5014554973288337,210.15458639016347,990.5911315659903,187.128,49.31274 +277,1.5607094914464252,196.63582783448533,1031.674709954466,188.231,47.63075 +278,1.5780719100901788,204.36109617554143,969.6154790568083,175.165,48.5366 +279,1.462355404138569,205.52826976961236,973.8990395981568,188.404,48.81386 +280,1.5713742641842252,219.18693337743758,973.1248619967882,176.469,50.27467 +281,1.5952530107032235,210.17788238492312,1140.2952821360864,204.013,49.20166 +282,1.4953350675173254,210.6552891302307,1029.2379578114267,195.145,49.37889 +283,1.6589426839683767,214.25642787422393,943.4188376038829,163.001,49.59871 +284,1.548449893513793,208.1908373145783,960.8199853679033,176.547,49.02505 +285,1.5454300679946056,210.15542793885578,1022.2860731609445,188.171,49.25949 +286,1.5006992610932648,211.68291563302716,1033.997271794892,195.417,49.49265 +287,1.4031516533647292,213.33286334304768,921.438446413576,185.047,49.8044 +288,1.4278850432220296,209.7079437807066,1006.7184802651564,198.998,49.34948 +289,1.4653246917713858,214.1551592883399,892.6259866903274,172.365,49.82415 +290,1.4796528125613722,207.075629325335,996.2730183728529,190.698,48.97643 +291,1.513943214232362,204.96053525188228,1057.9426141964982,198.368,48.68458 +292,1.5031501184167348,200.54026873344046,843.3126847682375,159.145,48.16953 +293,1.4993840353535344,208.89516341106824,991.7263909140294,187.576,49.16727 +294,1.498600877173489,211.1277991975509,1075.8119687878639,203.575,49.43027 +295,1.444444754954122,204.14014203549604,912.5340913690094,178.509,48.67015 +296,1.5144458333584354,216.22473689323223,965.4928449408392,180.979,50.00387 +297,1.3749150077920171,201.48003568694813,933.3570621763191,190.932,48.4346 +298,1.4721135556252394,207.95805686372637,1003.0864898384733,192.889,49.08977 +299,1.4844126834588776,215.24669805271705,945.9186155454905,180.536,49.92753 +300,1.4568841020262229,212.3943776926023,877.1486344059105,170.262,49.6292 +301,1.5194845852660002,210.51276591863135,988.2046008065466,184.684,49.33286 +302,1.5794138053574631,206.8889627735264,1164.8238536419262,210.27,48.83425 +303,1.4887284275385952,203.61491777514885,1005.593566479586,191.424,48.55461 +304,1.4461877399840617,206.2412207609325,1031.7626640012393,201.612,48.91788 +305,1.5277791044325062,206.8891614987287,1078.992411603712,200.667,48.89644 +306,1.5578753879010567,217.18150878281648,1042.3530725700577,190.49,50.0608 +307,1.568171435604285,200.5960078601297,974.567476109723,177.054,48.09914 +308,1.4566497371838287,210.59612930563696,954.2840624642611,185.262,49.41894 +309,1.4916921079025136,219.41087440102024,986.4224162863698,187.439,50.39914 +310,1.5206896158595107,208.83865335037882,1006.0570630365467,187.886,49.13485 +311,1.5003008340218036,212.24162032986578,1093.2246028437742,206.659,49.55841 +312,1.5823696998801666,211.54810040954584,1047.1637508828483,188.714,49.37746 +313,1.4424567745407217,209.5606162879219,1075.223157190687,210.596,49.31449 +314,1.5432984319208576,211.05804866397526,1070.911516765436,197.365,49.36775 +315,1.524538335449691,216.0015312912966,1023.1487325565689,190.645,49.96564 +316,1.4822150842214483,211.67920628466513,765.5524499317779,146.307,49.51473 +317,1.5191442868192524,191.75793587392607,862.1229274075665,161.153,47.08443 +318,1.5225884300470494,212.387646686334,907.0056975888056,169.198,49.54827 +319,1.464394897906053,216.85270494307065,961.2775894741859,185.728,50.13811 +320,1.6121460156201493,209.99040073411825,1032.6428941505014,183.023,49.15922 +321,1.5808437199249303,210.6551355809671,919.8124815447359,165.907,49.27499 +322,1.4945136274152855,209.82194631904957,890.8206678238748,168.985,49.28212 +323,1.4872083964329588,206.28903953001281,1073.9486205749947,204.624,48.87424 +324,1.5017305937092629,208.58602125426214,991.5054455501939,187.27,49.12804 +325,1.5111561009203085,213.0840709984489,920.0230217064495,172.793,49.6434 +326,1.4479206183689364,215.52339607701262,967.0130502236216,188.755,50.00445 +327,1.4291578953701565,212.21906354691552,1037.069584482615,204.832,49.64251 +328,1.4234642371371646,211.14740857731803,1058.84664739427,209.89,49.52394 +329,1.5187193225853133,213.15990741880697,887.65847215274,165.968,49.643 +330,1.5165520077498098,208.48690120264777,961.3767873020284,179.982,49.09845 +331,1.5435277804709864,217.39637255337115,1143.5496043680287,210.724,50.10327 +332,1.5657403935752954,209.4513290918491,1028.2176130421074,187.06,49.15229 +333,1.5383411435892123,210.60588974649733,1010.4611350329235,186.763,49.32086 +334,1.5411876565938092,209.70358676987033,941.6570450901071,173.757,49.21164 +335,1.479276346708207,213.1695518886687,1000.2394632959964,191.501,49.69232 +336,1.517010498843417,213.35039621101095,970.0126829681913,181.55,49.66726 +337,1.5378124081814903,211.89421301729197,901.2149530470817,166.622,49.47213 +338,1.5721152395593123,207.49627628121542,867.5923725451853,157.265,48.91467 +339,1.447747814520685,204.8947741966084,1004.1738654926027,196.03,48.75607 +340,1.486307474722688,208.31584016179704,1049.5761572385984,200.09,49.11484 +341,1.5169937680596353,211.58023935322143,1032.7862765862953,193.301,49.46081 +342,1.5308067334001119,199.4602332911128,970.8449398188432,180.233,48.00695 +343,1.3909894967755263,210.180299633989,999.6293622220524,202.341,49.44979 +344,1.5522766137394841,215.93817196360558,1016.2220038601894,186.314,49.92419 +345,1.4749819461301552,212.4258055627287,957.766589014055,183.851,49.61079 +346,1.4999303674278326,214.93068158480384,955.3782923208461,180.642,49.87182 +347,1.540283931962125,213.3870068434615,923.6166546295343,170.518,49.64307 +348,1.54869470641087,213.4928317268575,1008.1540709898126,185.219,49.64509 +349,1.4618838270986174,205.70396990712217,956.9869938477582,185.186,48.83528 +350,1.527129723164598,205.221874162993,1035.0009798895896,192.559,48.69979 +351,1.4379122235560804,208.97965707780577,1008.7718203659207,198.145,49.25158 +352,1.6222565072529458,203.11342172239816,1117.7611653192107,197.008,48.3355 +353,1.507492120699364,212.0558595900798,1009.4678330710564,190.007,49.52795 +354,1.5222363409097617,209.2315744781899,1116.3240145165967,208.289,49.17917 +355,1.477176442788764,215.66117286089624,948.4506303829581,181.819,49.98447 +356,1.4969874076799359,204.6401872385489,1036.4018976278016,196.308,48.66682 +357,1.5722562978463268,209.61321910089572,984.2491441036873,178.397,49.16339 +358,1.4442691518466035,208.85353673223847,872.5764592222778,170.711,49.22903 +359,1.4543561006640586,211.13393023054567,1105.5210558193116,214.928,49.48479 +360,1.6404605430614752,207.06707457464552,1114.8240125316352,194.546,48.78173 +361,1.5121876240200824,209.24376176626814,1034.6075561380921,194.195,49.19277 +362,1.503095824183998,210.5519786280611,967.3333469738665,182.555,49.35735 +363,1.4300945170100574,202.2818739873896,1052.259590662498,207.709,48.46521 +364,1.4861649954957143,207.01955550171917,1000.6804501378523,190.785,48.96196 +365,1.5397466126457813,210.90181971513198,947.5156050687226,174.985,49.35379 +366,1.5613146477824476,217.89592751490818,1039.5645324258282,189.605,50.13883 +367,1.5486831755672046,216.20979998406906,1090.407945312361,200.332,49.96 +368,1.5307049588556314,212.99173523537442,1041.5728458677202,193.375,49.60877 +369,1.445133116272648,213.50361578429795,891.5348823098757,174.326,49.77299 +370,1.4602700272437472,216.51782533904043,1067.2176130155358,206.722,50.10446 +371,1.574369914484293,208.64162808156354,907.8200926978586,164.346,49.04676 +372,1.4879781401464829,210.51015648620287,1006.3190973072008,191.649,49.37081 +373,1.556684023937835,208.19364787498276,1028.1627577442339,188.025,49.01544 +374,1.4505686425282698,212.25158639021905,1000.2486840148509,194.921,49.62021 +375,1.4625093330857741,208.82085318912573,985.2659402787823,190.584,49.20312 +376,1.4864922435331454,218.8567062113384,992.1243198153696,189.116,50.34189 +377,1.487805194491133,199.51274521194927,928.4634118416418,176.841,48.06411 +378,1.4868463287771823,206.5707366191168,1110.5366074702238,211.642,48.90804 +379,1.4917531112301892,216.44222124659126,1004.9363707240294,190.95,50.05695 +380,1.5018613918210788,209.5560177614375,1000.3037001297167,188.917,49.24198 +381,1.585939750152493,198.4019463662439,839.8763404449057,151.053,47.81442 +382,1.4728414584471945,210.07263940968784,1055.208353981121,202.821,49.33784 +383,1.5236543599013739,208.07877043559714,1081.523495459341,201.627,49.04179 +384,1.5515097701576366,210.92028915635782,1077.0126057989473,197.547,49.34165 +385,1.5136322631716064,207.19258486608558,1055.8253471200003,198.007,48.94932 +386,1.5055306561289339,210.00669923377242,1117.8362857751658,210.651,49.29045 +387,1.3774563935729958,202.91157914920214,905.0343658900504,184.828,48.60333 +388,1.4818369676215597,208.72695915116986,879.4139481781563,168.106,49.16869 +389,1.4507420983844586,205.69579160738215,1196.5939273660713,233.158,48.84769 +390,1.5900870973873573,207.94125291157675,970.9133367754561,174.213,48.9454 +391,1.4425398001405854,200.45934793387366,990.6495850873181,194.021,48.23163 +392,1.4501074008377988,209.57574734168995,871.9835320765218,169.974,49.307 +393,1.461901277084487,210.74313547651417,1029.296767822715,199.176,49.42981 +394,1.5208735979371006,203.3975550792148,1068.374858410845,199.503,48.49032 +395,1.4599894151699073,213.16428957303654,999.1489058705446,193.571,49.71527 +396,1.4311806250461026,208.41901147709405,1029.2290993651013,203.023,49.19361 +397,1.53064826902373,209.2950605127779,1122.5423112480842,208.415,49.17645 +398,1.550375538549328,208.71672941800122,1060.0097449709556,194.556,49.0846 +399,1.5966896742594094,219.80242300391512,901.7600095342002,161.207,50.31379 +400,1.4536733444421426,217.42459393366926,915.337515932989,178.029,50.21741 +401,1.4338729433354163,207.56902934650475,1033.6181367579668,203.543,49.08995 +402,1.4948801646694452,206.74718189797204,946.0630030620392,179.424,48.91925 +403,1.4643831467793647,205.6024518567053,1182.530630182085,228.478,48.82023 +404,1.5148755535592777,214.43379527202367,960.8738025504861,180.067,49.79582 +405,1.4907021661406976,206.43223009314823,787.8633332064638,149.798,48.887 +406,1.4310995361576067,200.94690660009803,820.0172977638714,161.763,48.30383 +407,1.440232051285095,216.62477839690814,996.4451065836516,195.439,50.14152 +408,1.5645281881160586,212.05491288450438,1062.6047060572128,193.45,49.45832 +409,1.4549551502347713,212.61975477096334,1058.2971387259477,205.671,49.65788 +410,1.4711676353883154,212.79770185308084,1045.2950206859866,201.122,49.65886 +411,1.4725223689434959,204.38615783864503,905.2697334481383,174.035,48.66587 +412,1.504325277096073,211.63504934394658,936.5979335633668,176.625,49.48264 +413,1.5013292670663048,212.85268518069242,1049.455257173224,198.263,49.62844 +414,1.5775940060061244,216.3849041609902,919.8966853794765,166.228,49.94463 +415,1.584179310476901,207.34511747722144,1083.5898249963543,195.079,48.88231 +416,1.5033993313025544,207.93925310771516,855.3842512180054,161.399,49.0498 +417,1.4914217744891827,202.20308382347196,1077.3839388674219,204.756,48.38278 +418,1.5758396328569624,206.01544505351654,1024.1576594015417,185.252,48.73534 +419,1.4637304238996014,206.85757634931605,1061.9595763092298,205.265,48.96981 +420,1.4980730869078294,218.69824020965507,961.4214017037664,181.987,50.30933 +421,1.5499668166237244,210.4870670647889,1063.6622407795387,195.273,49.29283 +422,1.4776421171335028,209.7522058122421,990.7600896560853,189.875,49.29438 +423,1.4581971442577106,213.31040767290875,1052.0199204462692,204.04,49.7345 +424,1.5296777347452784,217.68482826958822,1115.9475840989044,207.309,50.1536 +425,1.4812527655075431,211.8652615921709,950.111944859343,181.685,49.53765 +426,1.5011677658057276,207.00125275440567,1129.605624730884,213.426,48.94173 +427,1.518184907517148,208.90985435947786,1095.48948574851,204.892,49.14625 +428,1.3914692295102293,199.92948702567355,1063.667092417582,215.236,48.22828 +429,1.5030234224058683,211.2522315013038,893.888311734516,168.702,49.43945 +430,1.4988383144574817,220.0175470520067,1065.935969583669,201.678,50.4599 +431,1.5084853272571406,207.8047304494753,968.4341718504485,182.175,49.02779 +432,1.3935401055399992,211.2882021046058,1067.6386486785484,215.749,49.57685 +433,1.4212190336414692,206.6929155384659,907.4028294213263,180.128,49.00146 +434,1.451732440259562,210.2674081341095,833.247827600922,162.259,49.38632 +435,1.5307135219178534,213.7692753496684,979.911792957047,181.927,49.69922 +436,1.375964574500773,205.01790686998802,1098.6321547282157,224.585,48.85673 +437,1.3611002859635715,208.26238293903242,905.5999377226193,186.963,49.25975 +438,1.5402285150792419,206.62443023928267,1032.316321658849,190.593,48.85017 +439,1.5282730209584547,219.1811291986501,983.7646445059704,182.904,50.32742 +440,1.5188022184811307,202.57472260436282,1095.2330850634994,204.769,48.39461 +441,1.4619745575964318,207.19622917956127,998.2201764120631,193.154,49.01199 +442,1.4218557827330076,211.1438501446447,1077.397301834001,213.786,49.52548 +443,1.5146466911748646,202.0555121791064,990.0952564536273,185.569,48.33749 +444,1.5371014451489744,200.61005617017048,914.0187188945829,169.06,48.13766 +445,1.5438271802040762,211.3504771709864,1032.9577265701403,190.312,49.40129 +446,1.5929747403811656,212.342763892531,1052.6310408346544,188.57,49.45718 +447,1.6355372763609164,208.93097720954754,936.1693789105204,163.807,49.00675 +448,1.5213898199123568,212.1949551469445,957.8627686443033,178.812,49.52725 +449,1.4483801490222576,221.5492221104504,1026.1474776310417,200.241,50.69809 +450,1.4633857309298963,214.57965151139797,991.2683878750933,191.642,49.87588 +451,1.4873559276841084,212.81136435957248,896.6494185414648,170.828,49.64069 +452,1.5176313601958586,219.56477403865875,1037.2942129757848,194.071,50.38464 +453,1.4811893948108268,207.46447359785878,967.1728021323912,184.955,49.02054 +454,1.4873018660976112,208.11096733620866,1073.3049152859653,204.49,49.08948 +455,1.413339544839551,214.31000025206066,969.6250779885951,193.451,49.90585 +456,1.5188473892847687,204.36506044899568,1106.290485347955,206.831,48.60794 +457,1.5622817480085935,214.79795262875038,878.8196179365168,160.198,49.77993 +458,1.5414560535533026,212.02427543817856,871.5783002293946,160.801,49.48287 +459,1.4987851308114593,211.82267244516632,1044.486021905852,197.626,49.51132 +460,1.5088630683474427,211.4561570842369,1117.8889113802477,210.242,49.4562 +461,1.5613016954955339,195.66472523108627,974.6049156129855,177.758,47.51229 +462,1.4917674486228025,201.23950713636515,1007.534644658585,191.442,48.26695 +463,1.5624378886205252,204.56716070204706,1026.9120750951897,187.177,48.57978 +464,1.454208826135178,206.64043988160162,1118.191490733603,217.411,48.95555 +465,1.4728721089944088,212.31106244505312,850.302550715121,163.433,49.59997 +466,1.593115878069744,213.0132622315325,985.4000125879791,176.512,49.53503 +467,1.5287773254999948,202.95090285526095,1094.6024593402126,203.451,48.42763 +468,1.4838007612251785,208.9796524247578,960.6675508027613,183.419,49.19606 +469,1.4709674298627209,213.47335305814204,1020.4725480083379,196.37,49.73788 +470,1.5325874070701522,209.95602655257215,1041.0559848861744,193.066,49.25169 +471,1.4790575628037481,221.14865920704773,874.5655666855602,167.463,50.61406 +472,1.5630468661781491,216.83953725360485,931.8272955227525,169.786,50.015 +473,1.6114636805277986,214.9774514147266,1077.691120545833,191.08,49.74037 +474,1.501401109351294,208.91386920890662,983.2472938368313,185.747,49.16703 +475,1.519526529815976,217.0059313751695,1014.5051374289136,189.595,50.08785 +476,1.4544278188431357,214.1221353119249,1073.6920597115834,208.731,49.83365 +477,1.5327329764223894,213.72749453262492,978.3257495873753,181.417,49.6919 +478,1.453550541643471,205.65903284679612,869.1629401233713,169.062,48.83995 +479,1.5118474305772975,209.58275608386703,964.8009836376233,181.129,49.23302 +480,1.3522190587827465,203.76341265152914,1089.6672521122534,226.308,48.73539 +481,1.4667686999968435,211.8029583882842,888.7451656019795,171.463,49.54801 +482,1.5719352741179615,209.2798616997266,1023.2032317186828,185.491,49.12467 +483,1.6162889895940842,217.58220050294702,975.4747502970354,172.495,50.03484 +484,1.4703389907564821,208.37227494240236,1040.7918598781228,200.357,49.14078 +485,1.4947387932874565,208.02377927498264,1117.3745000567233,211.932,49.07022 +486,1.494289505411108,221.49794546518658,1007.2594051001738,191.098,50.63504 +487,1.5892659197112364,211.78659711355314,1077.0368458599557,193.344,49.39689 +488,1.4637694474241139,205.73120388706292,1036.776476401199,200.392,48.83625 +489,1.5494276621232659,210.20128251000233,897.0885331933042,164.744,49.26001 +490,1.4319870816855338,202.86431798906878,858.5902309648495,169.277,48.53267 +491,1.5054742679723545,215.62521960854014,998.9267370460153,188.25,49.94552 +492,1.3869792621331074,209.04619769832726,986.9760802693914,200.304,49.32105 +493,1.604799937511272,202.7143104354022,968.7457917920643,172.4,48.30879 +494,1.5852926340272282,212.24989170673274,1045.065045881314,188.025,49.45573 +495,1.505729221851194,211.6497176932863,1062.208762422283,200.145,49.48265 +496,1.4595256840553619,214.39029289454098,890.1597586491569,172.505,49.8586 +497,1.516418398641577,210.89370481643164,950.6965529446313,177.997,49.3812 +498,1.4988657362337177,206.88729923667225,914.2638938510604,172.978,48.93103 +499,1.4823192629193178,204.4800568231417,1025.8237971274598,196.036,48.66533 +500,1.442495725292382,213.55988511495525,905.4816745656944,177.345,49.78278 diff --git a/pyproject.toml b/pyproject.toml index 90eb451..e514ca7 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -48,6 +48,5 @@ line-length = 100 target-version = "py311" [tool.mypy] -python_version = "3.11" strict = true mypy_path = "src" \ No newline at end of file diff --git a/scripts/run_stochastic_study.py b/scripts/run_stochastic_study.py new file mode 100644 index 0000000..1f99275 --- /dev/null +++ b/scripts/run_stochastic_study.py @@ -0,0 +1,279 @@ +from pathlib import Path + +import matplotlib.pyplot as plt +import numpy as np +import pandas as pd + +from bodysimpy.analysis.stochastic import ( + StochasticStudyResult, + run_stochastic_study, +) +from bodysimpy.config.loader import load_config +from bodysimpy.modeling.crossmember import ( + build_crossmember_model, +) + + +def build_dataframe( + result: StochasticStudyResult, +) -> pd.DataFrame: + return pd.DataFrame( + { + "sample": [sample.sample_index for sample in result.samples], + "thickness_mm": [sample.thickness_m * 1000.0 for sample in result.samples], + "youngs_modulus_gpa": [sample.youngs_modulus_pa / 1e9 for sample in result.samples], + "tip_force_n": [sample.tip_force_n for sample in result.samples], + "max_stress_mpa": [sample.max_stress_pa / 1e6 for sample in result.samples], + "mode_1_frequency_hz": [sample.mode_1_frequency_hz for sample in result.samples], + } + ) + + +def save_histogram( + values: pd.Series, + *, + x_label: str, + title: str, + path: Path, +) -> None: + path.parent.mkdir( + parents=True, + exist_ok=True, + ) + + plt.figure(figsize=(8, 5)) + plt.hist(values, bins=30) + plt.xlabel(x_label) + plt.ylabel("Count") + plt.title(title) + plt.tight_layout() + plt.savefig(path, dpi=200) + plt.close() + + +def save_correlation_matrix( + dataframe: pd.DataFrame, + path: Path, +) -> None: + columns = [ + "thickness_mm", + "youngs_modulus_gpa", + "tip_force_n", + "max_stress_mpa", + "mode_1_frequency_hz", + ] + + values = dataframe[columns].to_numpy(dtype=float) + + correlation = np.corrcoef( + values, + rowvar=False, + ) + + path.parent.mkdir( + parents=True, + exist_ok=True, + ) + + plt.figure(figsize=(9, 7)) + + image = plt.imshow( + correlation, + vmin=-1.0, + vmax=1.0, + ) + + plt.colorbar( + image, + label="Pearson correlation", + ) + + labels = [ + "Thickness", + "E", + "Load", + "Stress", + "Mode 1", + ] + + plt.xticks( + range(len(labels)), + labels, + rotation=45, + ha="right", + ) + + plt.yticks( + range(len(labels)), + labels, + ) + + for row in range(len(labels)): + for column in range(len(labels)): + plt.text( + column, + row, + f"{correlation[row, column]:.2f}", + ha="center", + va="center", + ) + + plt.title("BodySimPy Monte Carlo Correlation Matrix") + + plt.tight_layout() + plt.savefig(path, dpi=200) + plt.close() + + +def save_sensitivity_plot( + sensitivities, + *, + title: str, + path: Path, +) -> None: + path.parent.mkdir( + parents=True, + exist_ok=True, + ) + + names = [point.parameter for point in sensitivities] + + values = [point.coefficient for point in sensitivities] + + plt.figure(figsize=(8, 5)) + plt.bar(names, values) + + plt.axhline( + 0.0, + linewidth=1.0, + ) + + plt.ylabel("Standardized regression coefficient") + + plt.title(title) + plt.tight_layout() + plt.savefig(path, dpi=200) + plt.close() + + +def main() -> None: + config = load_config("configs/stochastic_crossmember.yaml") + + if config.stochastic is None: + raise ValueError("Stochastic configuration is required.") + + model = build_crossmember_model(config) + + result = run_stochastic_study( + model, + config.stochastic, + ) + + dataframe = build_dataframe(result) + + validation_directory = Path("docs/validation") + + figure_directory = Path("docs/figures") + + validation_directory.mkdir( + parents=True, + exist_ok=True, + ) + + dataframe.to_csv( + validation_directory / "stochastic_samples.csv", + index=False, + ) + + correlation_columns = [ + "thickness_mm", + "youngs_modulus_gpa", + "tip_force_n", + "max_stress_mpa", + "mode_1_frequency_hz", + ] + + dataframe[correlation_columns].corr().to_csv( + validation_directory / "stochastic_correlation.csv" + ) + + save_histogram( + dataframe["max_stress_mpa"], + x_label="Maximum axial stress [MPa]", + title="Monte Carlo Stress Distribution", + path=(figure_directory / "stress_histogram.png"), + ) + + save_histogram( + dataframe["mode_1_frequency_hz"], + x_label="Mode-1 natural frequency [Hz]", + title="Monte Carlo Frequency Distribution", + path=(figure_directory / "frequency_histogram.png"), + ) + + save_correlation_matrix( + dataframe, + figure_directory / "stochastic_correlation_matrix.png", + ) + + save_sensitivity_plot( + result.stress_sensitivity, + title="Stress Sensitivity Ranking", + path=(figure_directory / "stress_sensitivity_ranking.png"), + ) + + save_sensitivity_plot( + result.frequency_sensitivity, + title="Mode-1 Frequency Sensitivity Ranking", + path=(figure_directory / "frequency_sensitivity_ranking.png"), + ) + + print() + print("BodySimPy Stochastic Engineering Study") + print("=" * 60) + + print() + print(f"Monte Carlo samples: {len(result.samples)}") + + print() + print("Maximum stress") + print("-" * 40) + print(f"Mean: {result.stress_summary.mean / 1e6:.4f} MPa") + print(f"Std. dev.: {result.stress_summary.standard_deviation / 1e6:.4f} MPa") + print(f"5th pct.: {result.stress_summary.percentile_5 / 1e6:.4f} MPa") + print(f"95th pct.: {result.stress_summary.percentile_95 / 1e6:.4f} MPa") + print(f"CoV: {result.stress_summary.coefficient_of_variation_percent:.3f} %") + print(f"P(stress > threshold): {result.stress_exceedance_probability_percent:.3f} %") + + print() + print("Mode-1 frequency") + print("-" * 40) + print(f"Mean: {result.frequency_summary.mean:.4f} Hz") + print(f"Std. dev.: {result.frequency_summary.standard_deviation:.4f} Hz") + print(f"5th pct.: {result.frequency_summary.percentile_5:.4f} Hz") + print(f"95th pct.: {result.frequency_summary.percentile_95:.4f} Hz") + print(f"CoV: {result.frequency_summary.coefficient_of_variation_percent:.3f} %") + + print() + print("Stress sensitivity") + print("-" * 40) + + for rank, item in enumerate( + result.stress_sensitivity, + start=1, + ): + print(f"{rank}. {item.parameter:<18} {item.coefficient:+.4f}") + + print() + print("Frequency sensitivity") + print("-" * 40) + + for rank, item in enumerate( + result.frequency_sensitivity, + start=1, + ): + print(f"{rank}. {item.parameter:<18} {item.coefficient:+.4f}") + + +if __name__ == "__main__": + main() diff --git a/src/bodysimpy/analysis/stochastic.py b/src/bodysimpy/analysis/stochastic.py new file mode 100644 index 0000000..4267e67 --- /dev/null +++ b/src/bodysimpy/analysis/stochastic.py @@ -0,0 +1,367 @@ +from __future__ import annotations + +from dataclasses import dataclass, replace + +import numpy as np +from numpy.typing import NDArray + +from bodysimpy.config.models import StochasticConfig +from bodysimpy.domain.structural_model import StructuralModel +from bodysimpy.solvers.calculix import CalculiXSolver + + +@dataclass(frozen=True, slots=True) +class StochasticInput: + """One randomized Monte Carlo input sample.""" + + sample_index: int + thickness_m: float + youngs_modulus_pa: float + tip_force_n: float + + +@dataclass(frozen=True, slots=True) +class DistributionSummary: + """Summary statistics for a sampled distribution.""" + + mean: float + standard_deviation: float + percentile_5: float + percentile_95: float + + +@dataclass(frozen=True, slots=True) +class SensitivityCoefficient: + """Standardized linear sensitivity coefficient.""" + + parameter: str + coefficient: float + + +@dataclass(frozen=True, slots=True) +class StochasticSample: + """Inputs and FEA outputs for one Monte Carlo realization.""" + + sample_index: int + thickness_m: float + youngs_modulus_pa: float + tip_force_n: float + max_stress_pa: float + mode_1_frequency_hz: float + + +@dataclass(frozen=True, slots=True) +class StochasticStudyResult: + """Aggregated result of a Monte Carlo FEA study.""" + + samples: tuple[StochasticSample, ...] + stress_summary: DistributionSummary + frequency_summary: DistributionSummary + stress_exceedance_probability_percent: float + stress_sensitivity: tuple[SensitivityCoefficient, ...] + frequency_sensitivity: tuple[SensitivityCoefficient, ...] + + +def _sample_positive_normal( + rng: np.random.Generator, + *, + mean: float, + standard_deviation: float, +) -> float: + """Draw from a normal distribution while enforcing positivity.""" + + while True: + value = float( + rng.normal( + loc=mean, + scale=standard_deviation, + ) + ) + + if value > 0.0: + return value + + +def generate_stochastic_inputs( + config: StochasticConfig, + *, + sample_count: int | None = None, +) -> tuple[StochasticInput, ...]: + """Generate reproducible randomized engineering inputs.""" + + count = config.samples if sample_count is None else sample_count + + if count <= 0: + raise ValueError("Stochastic sample count must be positive.") + + rng = np.random.default_rng(config.seed) + + samples: list[StochasticInput] = [] + + for sample_index in range( + 1, + count + 1, + ): + thickness_m = _sample_positive_normal( + rng, + mean=config.thickness_m.mean, + standard_deviation=(config.thickness_m.standard_deviation), + ) + + youngs_modulus_pa = _sample_positive_normal( + rng, + mean=config.youngs_modulus_pa.mean, + standard_deviation=(config.youngs_modulus_pa.standard_deviation), + ) + + tip_force_n = _sample_positive_normal( + rng, + mean=config.tip_force_n.mean, + standard_deviation=(config.tip_force_n.standard_deviation), + ) + + samples.append( + StochasticInput( + sample_index=sample_index, + thickness_m=thickness_m, + youngs_modulus_pa=youngs_modulus_pa, + tip_force_n=tip_force_n, + ) + ) + + return tuple(samples) + + +def summarize_distribution( + values: NDArray[np.float64], +) -> DistributionSummary: + """Calculate summary statistics for a sampled distribution.""" + + if values.ndim != 1: + raise ValueError("Distribution values must be one-dimensional.") + + if values.size == 0: + raise ValueError("Distribution values must not be empty.") + + if values.size > 1: + standard_deviation = float( + np.std( + values, + ddof=1, + ) + ) + else: + standard_deviation = 0.0 + + return DistributionSummary( + mean=float(np.mean(values)), + standard_deviation=standard_deviation, + percentile_5=float( + np.percentile( + values, + 5, + ) + ), + percentile_95=float( + np.percentile( + values, + 95, + ) + ), + ) + + +def calculate_sensitivity_ranking( + inputs: NDArray[np.float64], + output: NDArray[np.float64], + *, + parameter_names: tuple[str, ...], +) -> tuple[SensitivityCoefficient, ...]: + """Calculate standardized linear sensitivity coefficients.""" + + if inputs.ndim != 2: + raise ValueError("Sensitivity inputs must be a two-dimensional array.") + + if output.ndim != 1: + raise ValueError("Sensitivity output must be one-dimensional.") + + if inputs.shape[0] != output.shape[0]: + raise ValueError("Input and output sample counts must match.") + + if inputs.shape[1] != len(parameter_names): + raise ValueError("Parameter-name count must match input columns.") + + if inputs.shape[0] < 2: + raise ValueError("At least two samples are required for sensitivity analysis.") + + input_standard_deviation = np.std( + inputs, + axis=0, + ddof=1, + ) + + output_standard_deviation = float( + np.std( + output, + ddof=1, + ) + ) + + if np.any(input_standard_deviation == 0.0): + raise ValueError("Sensitivity inputs must vary.") + + if output_standard_deviation == 0.0: + raise ValueError("Sensitivity output must vary.") + + standardized_inputs = ( + inputs + - np.mean( + inputs, + axis=0, + ) + ) / input_standard_deviation + + standardized_output = (output - np.mean(output)) / output_standard_deviation + + coefficients, _, _, _ = np.linalg.lstsq( + standardized_inputs, + standardized_output, + rcond=None, + ) + + ranking = [ + SensitivityCoefficient( + parameter=name, + coefficient=float(coefficient), + ) + for name, coefficient in zip( + parameter_names, + coefficients, + strict=True, + ) + ] + + ranking.sort( + key=lambda item: abs(item.coefficient), + reverse=True, + ) + + return tuple(ranking) + + +def run_stochastic_study( + model: StructuralModel, + config: StochasticConfig, + *, + sample_count: int | None = None, +) -> StochasticStudyResult: + """Run Monte Carlo static and mode-1 FEA analyses.""" + + inputs = generate_stochastic_inputs( + config, + sample_count=sample_count, + ) + + solver = CalculiXSolver() + + samples: list[StochasticSample] = [] + + total = len(inputs) + + for input_sample in inputs: + section = replace( + model.section, + thickness_m=(input_sample.thickness_m), + ) + + material = replace( + model.material, + youngs_modulus_pa=(input_sample.youngs_modulus_pa), + ) + + sample_model = replace( + model, + name=(f"{model.name}_mc_{input_sample.sample_index:04d}"), + section=section, + material=material, + tip_force_n=(input_sample.tip_force_n), + ) + + static_result = solver.run(sample_model) + + modal_result = solver.run_modal( + sample_model, + modes=1, + ) + + if static_result.max_axial_stress_pa is None: + raise ValueError("Static FEA returned no axial stress.") + + samples.append( + StochasticSample( + sample_index=(input_sample.sample_index), + thickness_m=(input_sample.thickness_m), + youngs_modulus_pa=(input_sample.youngs_modulus_pa), + tip_force_n=(input_sample.tip_force_n), + max_stress_pa=(static_result.max_axial_stress_pa), + mode_1_frequency_hz=(modal_result.natural_frequencies_hz[0]), + ) + ) + + if ( + input_sample.sample_index == 1 + or input_sample.sample_index % 25 == 0 + or input_sample.sample_index == total + ): + print(f"Monte Carlo: {input_sample.sample_index}/{total}") + + stress_values = np.array( + [sample.max_stress_pa for sample in samples], + dtype=float, + ) + + frequency_values = np.array( + [sample.mode_1_frequency_hz for sample in samples], + dtype=float, + ) + + input_matrix = np.array( + [ + [ + sample.thickness_m, + sample.youngs_modulus_pa, + sample.tip_force_n, + ] + for sample in samples + ], + dtype=float, + ) + + exceedance_probability = float(np.mean(stress_values > config.stress_threshold_pa) * 100.0) + + parameter_names = ( + "thickness", + "youngs_modulus", + "tip_force", + ) + + return StochasticStudyResult( + samples=tuple(samples), + stress_summary=summarize_distribution(stress_values), + frequency_summary=summarize_distribution(frequency_values), + stress_exceedance_probability_percent=(exceedance_probability), + stress_sensitivity=( + calculate_sensitivity_ranking( + input_matrix, + stress_values, + parameter_names=parameter_names, + ) + ), + frequency_sensitivity=( + calculate_sensitivity_ranking( + input_matrix, + frequency_values, + parameter_names=parameter_names, + ) + ), + ) diff --git a/src/bodysimpy/config/models.py b/src/bodysimpy/config/models.py index 0a520ad..3406329 100644 --- a/src/bodysimpy/config/models.py +++ b/src/bodysimpy/config/models.py @@ -54,6 +54,20 @@ class AnalysisConfig(StrictModel): modal: ModalConfig +class NormalDistributionConfig(StrictModel): + mean: float + standard_deviation: float = Field(gt=0.0) + + +class StochasticConfig(StrictModel): + samples: int = Field(gt=0) + seed: int + stress_threshold_pa: float = Field(gt=0.0) + thickness_m: NormalDistributionConfig + youngs_modulus_pa: NormalDistributionConfig + tip_force_n: NormalDistributionConfig + + class SimulationConfig(StrictModel): project: ProjectConfig geometry: GeometryConfig @@ -61,3 +75,4 @@ class SimulationConfig(StrictModel): loading: LoadingConfig analysis: AnalysisConfig mesh: MeshConfig + stochastic: StochasticConfig | None = None diff --git a/tests/integration/test_stochastic_study.py b/tests/integration/test_stochastic_study.py new file mode 100644 index 0000000..f63b3a6 --- /dev/null +++ b/tests/integration/test_stochastic_study.py @@ -0,0 +1,38 @@ +import shutil + +import pytest + +from bodysimpy.analysis.stochastic import ( + run_stochastic_study, +) +from bodysimpy.config.loader import load_config +from bodysimpy.modeling.crossmember import ( + build_crossmember_model, +) + + +@pytest.mark.skipif( + shutil.which("ccx") is None, + reason="CalculiX is not installed.", +) +def test_stochastic_study_runs_real_solver() -> None: + config = load_config("configs/stochastic_crossmember.yaml") + + assert config.stochastic is not None + + model = build_crossmember_model(config) + + result = run_stochastic_study( + model, + config.stochastic, + sample_count=3, + ) + + assert len(result.samples) == 3 + + for sample in result.samples: + assert sample.max_stress_pa > 0.0 + assert sample.mode_1_frequency_hz > 0.0 + + assert result.stress_summary.mean > 0.0 + assert result.frequency_summary.mean > 0.0 diff --git a/tests/unit/test_config_loader.py b/tests/unit/test_config_loader.py index cce4a08..f725e47 100644 --- a/tests/unit/test_config_loader.py +++ b/tests/unit/test_config_loader.py @@ -1,5 +1,7 @@ from pathlib import Path +import pytest + from bodysimpy.config.loader import load_config @@ -18,3 +20,14 @@ def test_load_baseline_configuration() -> None: assert config.analysis.static is True assert config.analysis.modal.modes == 10 assert config.mesh.elements == 20 + + +def test_load_stochastic_configuration() -> None: + config = load_config("configs/stochastic_crossmember.yaml") + + assert config.stochastic is not None + assert config.stochastic.samples == 500 + assert config.stochastic.seed == 42 + assert config.stochastic.thickness_m.mean == pytest.approx(0.0015) + assert config.stochastic.thickness_m.standard_deviation == pytest.approx(0.00005) + assert config.stochastic.youngs_modulus_pa.mean == pytest.approx(210e9) diff --git a/tests/unit/test_stochastic.py b/tests/unit/test_stochastic.py new file mode 100644 index 0000000..0ae2e8e --- /dev/null +++ b/tests/unit/test_stochastic.py @@ -0,0 +1,135 @@ +import numpy as np +import pytest + +from bodysimpy.analysis.stochastic import ( + calculate_sensitivity_ranking, + generate_stochastic_inputs, + summarize_distribution, +) +from bodysimpy.config.models import StochasticConfig + + +def build_stochastic_config() -> StochasticConfig: + return StochasticConfig.model_validate( + { + "samples": 5, + "seed": 42, + "stress_threshold_pa": 350e6, + "thickness_m": { + "mean": 0.0015, + "standard_deviation": 0.00005, + }, + "youngs_modulus_pa": { + "mean": 210e9, + "standard_deviation": 5e9, + }, + "tip_force_n": { + "mean": 1000.0, + "standard_deviation": 80.0, + }, + } + ) + + +def test_distribution_summary() -> None: + values = np.array( + [ + 10.0, + 20.0, + 30.0, + 40.0, + 50.0, + ], + dtype=float, + ) + + summary = summarize_distribution(values) + + assert summary.mean == pytest.approx(30.0) + + assert summary.standard_deviation > 0.0 + + assert summary.percentile_5 < summary.mean + + assert summary.percentile_95 > summary.mean + + +def test_stochastic_inputs_are_reproducible() -> None: + config = build_stochastic_config() + + first = generate_stochastic_inputs(config) + + second = generate_stochastic_inputs(config) + + assert first == second + assert len(first) == 5 + + +def test_stochastic_inputs_are_physically_positive() -> None: + inputs = generate_stochastic_inputs(build_stochastic_config()) + + for sample in inputs: + assert sample.thickness_m > 0.0 + + assert sample.youngs_modulus_pa > 0.0 + + assert sample.tip_force_n > 0.0 + + +def test_sample_count_override() -> None: + inputs = generate_stochastic_inputs( + build_stochastic_config(), + sample_count=3, + ) + + assert len(inputs) == 3 + + assert [sample.sample_index for sample in inputs] == [ + 1, + 2, + 3, + ] + + +def test_sensitivity_ranking() -> None: + inputs = np.array( + [ + [1.0, 10.0, 100.0], + [2.0, 8.0, 110.0], + [3.0, 6.0, 90.0], + [4.0, 4.0, 120.0], + [5.0, 2.0, 80.0], + ], + dtype=float, + ) + + output = np.array( + [ + 2.0, + 4.0, + 6.0, + 8.0, + 10.0, + ], + dtype=float, + ) + + ranking = calculate_sensitivity_ranking( + inputs, + output, + parameter_names=( + "parameter_a", + "parameter_b", + "parameter_c", + ), + ) + + assert len(ranking) == 3 + + assert {item.parameter for item in ranking} == { + "parameter_a", + "parameter_b", + "parameter_c", + } + + assert abs(ranking[0].coefficient) >= abs(ranking[1].coefficient) >= abs(ranking[2].coefficient)