* Objective: convex quadratic * Constraints: nonconvex nonlinear *Knowing a list of points $(x_i, y_i)$, describing the center line of *a wood piece, we want to fit the best polynomial that minimizes the *sum of the squares of errors (least square problem). This part is *simple. But we have to deal with three sets of conditions : * *1) the polynomial must go througn the first point $(x_0, y_0)$ of * the list (easy...); *2) the initial slope of the curve $m_0$ is also known (easy too...); * *3) the radius of curvature at every point must exceed a given radius * $R$. In fact, I modified this last condition. I prefer choosing a * condition about the square of the radius to avoid square root and * absolute value functions. If this had not been satisfying, I would * have tried the exact condition. * *From: * *Francois Grondin, Ph.D. *Forintek Canada Corp. *319, rue Franquet, Quebec, QC *G1P 4R4 *Tel.: (418) 659-2647 *Fax : (418) 659-2922 *E-Mail : francois.grondin@qc.forintek.ca set j /0*195/; set i /0*4/; parameter R; R=2500; parameter x[j]; parameter y[j]; x['0']= 0.000000; y['0']= 3.556000; x['1']= 2.540000; y['1']= 3.563460; x['2']= 5.080000; y['2']= 3.570910; x['3']= 7.620010; y['3']= 3.578370; x['4']= 10.161200; y['4']= 3.593710; x['5']= 12.701200; y['5']= 3.602130; x['6']= 15.241200; y['6']= 3.610550; x['7']= 17.781200; y['7']= 3.618980; x['8']= 20.321200; y['8']= 3.627400; x['9']= 22.861200; y['9']= 3.635820; x['10']= 25.401200; y['10']= 3.642300; x['11']= 27.941200; y['11']= 3.650730; x['12']= 30.481200; y['12']= 3.665300; x['13']= 33.021200; y['13']= 3.673810; x['14']= 35.561200; y['14']= 3.682330; x['15']= 38.101200; y['15']= 3.690840; x['16']= 40.641100; y['16']= 3.506370; x['17']= 43.181100; y['17']= 3.514980; x['18']= 45.723100; y['18']= 2.774810; x['19']= 48.263100; y['19']= 2.783760; x['20']= 50.803100; y['20']= 2.792720; x['21']= 53.343100; y['21']= 2.801680; x['22']= 55.883900; y['22']= 2.785290; x['23']= 58.423900; y['23']= 2.794150; x['24']= 60.965000; y['24']= 2.476020; x['25']= 63.505000; y['25']= 2.486170; x['26']= 66.045000; y['26']= 2.483690; x['27']= 68.585000; y['27']= 2.493770; x['28']= 71.125000; y['28']= 2.276350; x['29']= 73.665000; y['29']= 2.285210; x['30']= 76.205000; y['30']= 2.267440; x['31']= 78.745000; y['31']= 2.276050; x['32']= 81.285000; y['32']= 2.284650; x['33']= 83.825000; y['33']= 2.293250; x['34']= 86.365000; y['34']= 2.490700; x['35']= 88.905000; y['35']= 2.499210; x['36']= 91.445000; y['36']= 2.507730; x['37']= 93.985000; y['37']= 2.516240; x['38']= 96.525000; y['38']= 2.533810; x['39']= 99.065000; y['39']= 2.542410; x['40']= 101.605000; y['40']= 2.748550; x['41']= 104.145000; y['41']= 2.757160; x['42']= 106.685000; y['42']= 2.765760; x['43']= 109.225000; y['43']= 2.774370; x['44']= 111.766000; y['44']= 2.607580; x['45']= 114.306000; y['45']= 2.616100; x['46']= 116.847000; y['46']= 2.260080; x['47']= 119.387000; y['47']= 2.268600; x['48']= 121.927000; y['48']= 2.286950; x['49']= 124.467000; y['49']= 2.295560; x['50']= 127.008000; y['50']= 1.478080; x['51']= 129.548000; y['51']= 1.485400; x['52']= 132.090000; y['52']= 1.110100; x['53']= 134.630000; y['53']= 1.118800; x['54']= 137.170000; y['54']= 1.127490; x['55']= 139.710000; y['55']= 1.136190; x['56']= 142.251000; y['56']= 0.978377; x['57']= 144.791000; y['57']= 0.987075; x['58']= 147.330000; y['58']= 0.822030; x['59']= 149.870000; y['59']= 0.829676; x['60']= 152.410000; y['60']= 0.837322; x['61']= 154.950000; y['61']= 0.844885; x['62']= 157.490000; y['62']= 0.852781; x['63']= 160.030000; y['63']= 0.860510; x['64']= 162.571000; y['64']= 0.538377; x['65']= 165.111000; y['65']= 0.546996; x['66']= 167.651000; y['66']= 0.600516; x['67']= 170.191000; y['67']= 0.609224; x['68']= 172.731000; y['68']= 0.789128; x['69']= 175.271000; y['69']= 0.798084; x['70']= 177.811000; y['70']= 1.009960; x['71']= 180.351000; y['71']= 1.018830; x['72']= 182.891000; y['72']= 1.037120; x['73']= 185.431000; y['73']= 1.045910; x['74']= 187.972000; y['74']= 0.720837; x['75']= 190.512000; y['75']= 0.699625; x['76']= 193.052000; y['76']= 0.708414; x['77']= 195.592000; y['77']= 0.717203; x['78']= 198.133000; y['78']= 0.609086; x['79']= 200.673000; y['79']= 0.617627; x['80']= 203.213000; y['80']= 0.199369; x['81']= 205.753000; y['81']= 0.207992; x['82']= 208.293000; y['82']= 0.225288; x['83']= 210.833000; y['83']= 0.233831; x['84']= 213.374000; y['84']= -0.172432; x['85']= 215.914000; y['85']= -0.163889; x['86']= 218.456000; y['86']= -0.348890; x['87']= 220.996000; y['87']= -0.341233; x['88']= 223.536000; y['88']= -0.333576; x['89']= 226.076000; y['89']= -0.325920; x['90']= 228.616000; y['90']= -0.522016; x['91']= 231.156000; y['91']= -0.514426; x['92']= 233.696000; y['92']= -0.506837; x['93']= 236.236000; y['93']= -0.499248; x['94']= 238.776000; y['94']= -0.491658; x['95']= 241.316000; y['95']= -0.484069; x['96']= 243.856000; y['96']= -0.476480; x['97']= 246.396000; y['97']= -0.468890; x['98']= 248.936000; y['98']= -0.462453; x['99']= 251.476000; y['99']= -0.455008; x['100']= 254.016000; y['100']= -0.447563; x['101']= 256.556000; y['101']= -0.440118; x['102']= 259.095000; y['102']= 0.055237; x['103']= 261.635000; y['103']= 0.063873; x['104']= 264.176000; y['104']= 0.062826; x['105']= 266.716000; y['105']= 0.070575; x['106']= 269.256000; y['106']= -0.275144; x['107']= 271.797000; y['107']= -0.275842; x['108']= 274.337000; y['108']= -0.268089; x['109']= 276.877000; y['109']= -0.260335; x['110']= 279.417000; y['110']= -0.118406; x['111']= 281.957000; y['111']= -0.110665; x['112']= 284.497000; y['112']= -0.102923; x['113']= 287.037000; y['113']= -0.095181; x['114']= 289.577000; y['114']= 0.103362; x['115']= 292.117000; y['115']= 0.111028; x['116']= 294.657000; y['116']= 0.131239; x['117']= 297.197000; y['117']= 0.138981; x['118']= 299.736000; y['118']= 0.272430; x['119']= 302.276000; y['119']= 0.280024; x['120']= 304.815000; y['120']= 0.421614; x['121']= 307.355000; y['121']= 0.429144; x['122']= 309.895000; y['122']= 0.436674; x['123']= 312.435000; y['123']= 0.444203; x['124']= 314.976000; y['124']= 0.484027; x['125']= 317.516000; y['125']= 0.491804; x['126']= 320.056000; y['126']= 0.515482; x['127']= 322.598000; y['127']= 0.571256; x['128']= 325.138000; y['128']= 0.580224; x['129']= 327.678000; y['129']= 0.589193; x['130']= 330.215000; y['130']= 1.111820; x['131']= 332.755000; y['131']= 1.120710; x['132']= 335.295000; y['132']= 1.345510; x['133']= 337.835000; y['133']= 1.354480; x['134']= 340.375000; y['134']= 1.383820; x['135']= 342.915000; y['135']= 1.392870; x['136']= 345.454000; y['136']= 1.772980; x['137']= 347.994000; y['137']= 1.782020; x['138']= 350.534000; y['138']= 2.006330; x['139']= 353.074000; y['139']= 2.015440; x['140']= 355.615000; y['140']= 2.412930; x['141']= 358.155000; y['141']= 2.420930; x['142']= 360.695000; y['142']= 2.390290; x['143']= 363.235000; y['143']= 2.398230; x['144']= 365.773000; y['144']= 2.891230; x['145']= 368.313000; y['145']= 2.899110; x['146']= 370.853000; y['146']= 3.337140; x['147']= 373.393000; y['147']= 3.345090; x['148']= 375.933000; y['148']= 3.404410; x['149']= 378.473000; y['149']= 3.412480; x['150']= 381.012000; y['150']= 3.702710; x['151']= 383.552000; y['151']= 3.710940; x['152']= 386.090000; y['152']= 4.153530; x['153']= 388.632000; y['153']= 4.127350; x['154']= 391.170000; y['154']= 4.560110; x['155']= 393.710000; y['155']= 4.569230; x['156']= 396.250000; y['156']= 4.772520; x['157']= 398.790000; y['157']= 4.781630; x['158']= 401.329000; y['158']= 4.891490; x['159']= 403.869000; y['159']= 4.899350; x['160']= 406.410000; y['160']= 4.964780; x['161']= 408.950000; y['161']= 4.972720; x['162']= 411.490000; y['162']= 5.167670; x['163']= 414.030000; y['163']= 5.175530; x['164']= 416.570000; y['164']= 5.385770; x['165']= 419.110000; y['165']= 5.393640; x['166']= 421.650000; y['166']= 5.401510; x['167']= 424.190000; y['167']= 5.409370; x['168']= 426.729000; y['168']= 5.283580; x['169']= 429.269000; y['169']= 5.291360; x['170']= 431.810000; y['170']= 5.142030; x['171']= 434.350000; y['171']= 5.148650; x['172']= 436.891000; y['172']= 4.800670; x['173']= 439.431000; y['173']= 4.807290; x['174']= 441.972000; y['174']= 4.861690; x['175']= 444.512000; y['175']= 4.869470; x['176']= 447.052000; y['176']= 4.561590; x['177']= 449.592000; y['177']= 4.569370; x['178']= 452.134000; y['178']= 3.788240; x['179']= 454.674000; y['179']= 3.796020; x['180']= 457.214000; y['180']= 3.803800; x['181']= 459.754000; y['181']= 3.811580; x['182']= 462.295000; y['182']= 3.665340; x['183']= 464.835000; y['183']= 3.673060; x['184']= 467.375000; y['184']= 3.676650; x['185']= 469.915000; y['185']= 3.684500; x['186']= 472.455000; y['186']= 3.690400; x['187']= 474.995000; y['187']= 3.698320; x['188']= 477.535000; y['188']= 3.706240; x['189']= 480.075000; y['189']= 3.714160; x['190']= 482.616000; y['190']= 3.501540; x['191']= 485.156000; y['191']= 3.509630; x['192']= 487.696000; y['192']= 3.517720; x['193']= 490.236000; y['193']= 3.525810; x['194']= 492.776000; y['194']= 3.533900; x['195']= 495.316000; y['195']= 3.541990; **************************************************** parameter m_0; m_0=(y['15']-y['0'])/(x['15']-x['0']); variable a[i], p[j], pprime[j], pprime2[j], z; equation first[j], second[j], third[j], curvature, defobj; first[j].. p[j]=e=a['0']+ a['1']*x[j]+ a['2']*x[j]*x[j]+ a['3']*x[j]*x[j]*x[j]+ a['4']*x[j]*x[j]*x[j]*x[j]; second[j].. pprime[j]=e=1*a['1']+ 2*a['2']*x[j]+ 3*a['3']*x[j]*x[j]+ 4*a['4']*x[j]*x[j]*x[j]; third[j].. pprime2[j]=e= 2*1*a['2']+ 3*2*a['3']*x[j]+ 4*3*a['4']*x[j]*x[j]; curvature[j].. sqr(R) * sqr(pprime2[j]) =l= ( 1 + sqr(pprime[j]))* ( 1 + sqr(pprime[j]))* ( 1 + sqr(pprime[j])); defobj.. z=e= sum(j, sqr(p[j] - y[j]) ); p.fx['0'] = y['0']; pprime.fx['0'] = m_0; model sawpath /all/; solve sawpath using dnlp minimaze z; file dem /output/; put dem; Put "z = ", z.l:20:10; Put /; Put "points p prime prime2 "; Put /; Put " ------------------------------------------------ "; Put /; Loop(j,put j.tl:3 Put p.l[j]:11:5 Put pprime.l[j]:11:5 Put pprime2.l[j]:11:5 Put /;); Put "points a[] "; Put /; Loop(i,put i.tl:3 Put a.l[i]:11:5 Put /;);