求参数方程的系数的问题(拟合)

laihai718

首先已知一系列点x，y的值，x，y都是t的函数（摆线），x，y的方程形式可写出如下：
摆线参数方程
x=a*sin(b+t)-c*sin(b+g*t);
y=a*cos(b+t)-c*cos(b+g*t);
现在需要求出系数a、b、c、g。其中，t是摆线的摆角。

已知数据点坐标如下：
           x                   y

37.856908	7.364225
37.722124	7.385029
37.58687	7.412964
37.329124	7.463645
36.845487	7.608982
36.824192	7.615679
36.393619	7.791825
36.302059	7.830047
36.058921	7.955129
35.8536	8.062848
35.699383	8.16122
35.458248	8.32491
35.384997	8.373934
35.197788	8.523615
35.028576	8.668206
34.99186	8.698158
34.96201	8.723733
34.82907	8.855671
34.749556	8.938238
34.621844	9.076793
34.570672	9.131525
34.382182	9.367993
34.296894	9.473604
34.18626	9.624569
34.16005	9.663583
34.143829	9.690129
34.072003	9.807464
33.963101	9.981147
33.889993	10.093898
33.846489	10.173587
33.776284	10.310817

原来用1stopt软件可以得到解，但现在不行了。不知道为什么。
下面是我原来的方法，从已知摆线上取10个点，得到的数据再代入下面的1stopt程序，就可以得到对应的系数。
SharedModel;
Variable t,x,y;
Function
x=a*sin(b+c+t)-e*cos(f)*sin(b+c+g*t);
y=a*cos(b+c+t)-e*cos(f)*cos(b+c+g*t);
Data;
0        13.7836        12.9607
0.1        13.9968        12.9333
0.2        14.4546        13.0332
0.3        15.1899        13.1918
0.4        16.2206        13.3350
0.5        17.5497        13.3855
0.6        19.1641        13.2660
0.7        21.0355        12.9017
0.8        23.1201        12.2234
0.9        25.3607        11.1702
1        27.6875        9.6919

刚刚搜了一下论坛，发现有个帖子跟我的问题很类似，只是他的是一个方程，我的是两个方程，http://forum.simwe.com/viewthread.php?tid=872573&highlight=%B7%BD%B3%CC%CF%B5%CA%FD

真的是不知道该如何办了，还希望大家帮帮我，谢谢。

messenger

参数方程拟合1stopt要用SharedModel模式。

具体参考，http://forum.simwe.com/thread-862684-1-1.html，http://forum.simwe.com/thread-810618-1-1.html。

不过估计需要1stopt 2.5版本。你这2个参数方程比较简单，应该可以合成1个方程，从而不用SharedModel模式也能解出。

laihai718

这两个参数方程好像不好合并，用1stopt来解得话，不稳定，能帮我指点一下matlab的方法吗

messenger

matlab的方法就是最小二乘法，最原始的方法也是适用范围最广的方法。

这两个参数方程好像不好合并，用1stopt来解得话，不稳定，能帮我指点一下matlab的方法吗
laihai718 发表于 2010-5-4 18:28

laihai718

苦恼~~，谢谢版主哈

shamohu

你的两组参数方程都是周期函数，不加约束，会有无数组解，也就是楼主说的“不稳定吧”。

第二组方程：
x=a*sin(b+c+t)-e*cos(f)*sin(b+c+g*t);
y=a*cos(b+c+t)-e*cos(f)*cos(b+c+g*t);

b+c完全可简化成一个啊！

lengyunfeng

LZ，如果我是你的话，我会再去找些数据，数据的量以达到X、Y实现数个周期为准，这样对应的你能找到周期T，再用(2*pi,2*pi/g)的最小公倍数=T来确定g；三角函数内的b是用来移动对称轴的，你找到min(y)，对应的找到t的值，按照摆线的理论形式，一般来说到达最低点y时，y'是=0的，你就可以依此找到b。b、g出来了，三角函数外的a、c就很好办了

lin2009

1# laihai718
Variable x,y;
ParVariable t[0,2*pi];
SharedModel;
Function x=a*sin(b+t)-c*sin(b+g*t);
         y=a*cos(b+t)-c*cos(b+g*t);
data;
37.856908 7.364225
37.722124 7.385029
...
迭代数: 23
计算用时(时:分:秒:微秒): 00:04:07:844
优化算法: 通用全局优化法(UGO1)
计算结束原因: 用户中止
均方差(RMSE): 0.00238806227306983
残差平方和(RSS): 0.000353576168043685
相关系数(R): 0.999999995699039
相关系数之平方(R^2): 0.999999991398078
决定系数(DC): 0.999999991398077
F统计(F-Statistic): -30048.2359821257
参数           最佳估算
-------------------- -------------
a           -5.27183652550824
b           -1.88424680219221
c           40.5319682246089
g           -0.00126258915569132
====== 输出结果 =====
文件：数据文件 - 1
No     目标x                   计算x                      目标y                计算y
1    37.856908     37.8571877000443      7.364225      7.36627214998034
2    37.722124     37.7223567050368      7.385029      7.38645800634685
3    37.58687      37.5863867279656      7.412964      7.41041515053898
4    37.329124     37.3296186198965      7.463645      7.46569526178463
5    36.845487     36.8448907101895      7.608982      7.60724965644107
6    36.824192     36.8238010026932      7.615679      7.61455833358351
7    36.393619     36.3914110301499      7.791825      7.78691937726332
8    36.302059     36.3013182691603      7.830047      7.82847970666481
9    36.058921     36.0570613121731      7.955129      7.95166497258738
10   35.8536      35.8552878658093      8.062848      8.06568924326164
11   35.699383    35.6997172534579      8.16122       8.16174123624509
12   35.458248    35.4587910473693      8.32491       8.32566322570448
13   35.384997    35.3878226024086      8.373934      8.37772044360328
14   35.197788    35.1991084630843      8.523615      8.52522864669617
15   35.028576    35.0298834037384      8.668206      8.66967481788107
16   34.99186     34.9946893044085      8.698158      8.70128069301
17   34.96201     34.9656970206221      8.723733      8.7277429520819
18   34.82907     34.8301726489544      8.855671      8.85678926526957
19   34.749556    34.74943418302        8.938238      8.93811962636532
20   34.621844    34.6204889763345      9.076793      9.07556580970654
21   34.570672    34.5705844680708      9.131525      9.13144792279341
22   34.382182    34.3771387754698      9.367993      9.36403346238077
23   34.296894    34.2946244331034      9.473604      9.4719131349948
24   34.18626     34.1854045258618      9.624569      9.62397688458574
25   34.16005     34.1588198080556      9.663583      9.66274735533961
26   34.143829    34.1414245118621      9.690129      9.6885159417782
27   34.072003    34.0662146683838      9.807464      9.80379238124653
28   33.963101    33.9596359959749      9.981147      9.97913040596496
29   33.889993    33.8936205873264      10.093898     10.0958896875907
30   33.846489    33.8506350712725      10.173587     10.1757730842536
31   33.776284    33.7807653272918      10.310817     10.3130181055728

laihai718

你的两组参数方程都是周期函数，不加约束，会有无数组解，也就是楼主说的“不稳定吧”。

第二组方程：
x=a*sin(b+c+t)-e*cos(f)*sin(b+c+g*t);
y=a*cos(b+c+t)-e*cos(f)*cos(b+c+g*t);

b+c完全可简化成一个啊 ...
shamohu 发表于 2010-5-4 22:03

谢谢哈，这样也可以

laihai718

LZ，如果我是你的话，我会再去找些数据，数据的量以达到X、Y实现数个周期为准，这样对应的你能找到周期T，再用max(2*pi,2*pi/g)=T来确定g；三角函数内的b是用来移动对称轴的，你找到min(y)，对应的找到t的值，按照摆 ...
lengyunfeng 发表于 2010-5-4 22:28

点数有限，不能取到一个周期的

laihai718

1# laihai718
Variable x,y;
ParVariable t[0,2*pi];
SharedModel;
Function x=a*sin(b+t)-c*sin(b+g*t);
         y=a*cos(b+t)-c*cos(b+g*t);
data;
37.856908 7.364225
37.722124 7.385029
...
迭代数 ...
lin2009 发表于 2010-5-4 22:36

您用的1stopt是哪个版本的？我的1.0的版本，运行不了

我的想法是用已知的x  y的数据，推出t的数据，然后根据这3组数据在1stopt里实现，之前就可以，现在确不行了