提高超越方程解的精度

Question

我有一个特定的运动学作为更复杂机器的一部分，并需要计算一些物理参数，这些参数使用我拥有的仪器很难（更像是不可能）以适当的精度进行测量。

[运动学]

乍一看，这是一个简单的1自由度臂（黑色），可以围绕x轴旋转。它有一个重量，强制它一直向上移动，直到它撞到机械端点（角度a0）或某个半径为r0的管子（蓝色）。臂旋转中心在y0处。管子可以移动到任何y(t)高度。

[用途]

这用于测量管子的半径以进一步处理。可以通过基本的角度测量计算出半径，从而得到图像底部的方程式。常数a0、y0、z0很难测量（它在复杂的机械内部），因此距离的测量精度最小为0.1毫米，角度为0.1度，甚至这也是值得质疑的。

[校准]

因此，我决定尝试从机器本身完成的一组测量中计算这些参数（自动校准）。因此，我有一个已知半径r0的校准管。所有绿色参数都可以处理为常数。现在，我将管子沿着y轴定位，以尽可能覆盖臂的角度。遗憾的是，范围仅约为20度（对于当前的机器设置），记住预设y(t)的测量a(t)...作为n点数据集。这给我提供了一个由n个超越方程组成的系统。从这里，我尝试/猜测“所有”可能性a0、y0、z0，并记住最佳解决方案（最接近r0）。

[a0、y0、z0的逼近]

逼近基于我的这个类：

//---------------------------------------------------------------------------
class approx
    {
public:
    double a,aa,a0,a1,da,*e,e0;
    int i,n;
    bool done,stop;

    approx()            { a=0.0; aa=0.0; a0=0.0; a1=1.0; da=0.1; e=NULL; e0=NULL; i=0; n=5; done=true; }
    approx(approx& a)   { *this=a; }
    ~approx()           {}
    approx* operator = (const approx *a) { *this=*a; return this; }
    //approx* operator = (const approx &a) { ...copy... return this; }

    void init(double _a0,double _a1,double _da,int _n,double *_e)
        {
        if (_a0<=_a1) { a0=_a0; a1=_a1; }
        else          { a0=_a1; a1=_a0; }
        da=fabs(_da);
        n =_n ;
        e =_e ;
        e0=-1.0;
        i=0; a=a0; aa=a0;
        done=false; stop=false;
        }
    void step()
        {
        if ((e0<0.0)||(e0>*e)) { e0=*e; aa=a; }         // better solution
        if (stop)                                       // increase accuracy
            {
            i++; if (i>=n) { done=true; a=aa; return; } // final solution
            a0=aa-fabs(da);
            a1=aa+fabs(da);
            a=a0; da*=0.1;
            a0+=da; a1-=da;
            stop=false;
            }
        else{
            a+=da; if (a>a1) { a=a1; stop=true; }       // next point
            }
        }
    };
//---------------------------------------------------------------------------

它通过一些初始步骤搜索单变量的全范围，然后找到最小偏差点。之后改变范围和步长以接近该点的区域，并递归提高精度。

解决方案本身如下图所示：

// (global) input data
#define _irc_calib_n 100
#define _irc_approx_n 5
int    irc_calib_ix; // number of measured points
double irc_calib_y[_irc_calib_n]; // y(t)
double irc_calib_a[_irc_calib_n]; // a(t)
double irc_calib_r; // calibration tube radius + arm radius

// approximation
int ix=0;
double e,a,deg=M_PI/180.0;
approx aa,ay,az;
//           min       max       step     recursions    ErrorOfSolutionVariable
for (aa.init(-90.0*deg,+90.0*deg,10.0*deg,_irc_approx_n,&e);!aa.done;aa.step())
for (ay.init(  0.0    ,200.0    ,10.0    ,_irc_approx_n,&e);!ay.done;ay.step())
for (az.init( 50.0    ,400.0    ,10.0    ,_irc_approx_n,&e);!az.done;az.step())
    {
    for (e=0.0,ix=0;ix<_irc_calib_n;ix++) // test all measured points (e is cumulative error)
        {
        a=irc_calib_a[ix]+aa.a;
        if (a> pi) a-=pi2;
        if (a<-pi) a+=pi2;
        if (fabs(a)>0.5*pi) { e=100.0; break; } // ignore too far angles
        e+=fabs(+(cos(a)*(irc_calib_y[ix]-ay.a))
                -(sin(a)*(az.a))
                -(irc_calib_r));
        }
    }
// here aa.a,ay.a,az.a holds the result

这导致解决方案接近测量值，但在模拟中结果仍然不够准确。这取决于点数和角度范围，误差在0.1毫米到0.5毫米之间。如果我正确测量z0并忽略其近似，则精度显着提高，使得y0没有误差（在模拟中），而a0的误差约为0.3度。

Q1：如何进一步提高解决方案的准确性？

我无法增加角度范围。点数最好在100左右，点数越多，准确性越好，但超过150时结果不稳定（对于某些半径完全偏离）。我不知道原因。递归次数大于6没有太大影响。

是否可以根据距离0度的角度距离加权偏差？但遗憾的是，a（t）范围不一定包括0度。

所需精度为0.01毫米，对于y0，z0和a0为0.01度

Q2：我有遗漏的地方吗？

例如，错误嵌套的近似或一些数学简化或不同的方法。

[注]

角度必须以a（t）+ a0的形式表示，因为它是通过带有SW重置的IRC（16000个步骤/回合）测量的。当在a0位置时重置。我不计算振动和校准管偏心，它们已经被照顾好了，我的第一个目标是使这个模拟工作。管y（t）可以任意定位，a（t）测量可以随意进行。

现在，校准过程沿着y轴扫描点（从a0向下移动）。使用6次递归的计算需要大约35秒（请耐心等待）。5次递归需要大约22秒。

[编辑1]以下是模拟的执行方式

approx aa; double e;
for (aa.init(-90.0*deg,+90.0*deg,10.0*deg,6,&e);!aa.done;aa.step())
 e=fabs(+(cos(aa.a)*(y(t)-y0))
        -(sin(aa.a)*(z0))
        -(irc_calib_r));
if (aa.a<a0) aa.a=a0;

[编辑2] 一些数值

刚刚意识到模拟代码中只有4次递归，为了匹配输入IRC的精度，必须进行6次递归。在更改后（也包括之前的编辑），以下是一些结果。

                | a0[deg]| y0[mm] | z0[mm] | 
    simulated   | -7.4510|191.2590|225.9000|
    z0 known    | -7.4441|191.1433|225.9000|
    z0 unknown  | -7.6340|191.8074|225.4971|

因此，使用已知的z0测量的准确度几乎在期望范围内，但是使用未知的z0时，误差仍然比所需大约10倍。增加模拟精度对于递归次数大于6没有影响，并且也没有意义，因为实际输入数据也不会更加准确。

这里是使用上述模拟设置进行测试的模拟/测量点：

 ix   a [deg]    y [mm]
  0   -0.2475 +105.7231 
  1   -0.4500 +104.9231 
  2   -0.6525 +104.1231 
  3   -0.8550 +103.3231 
  4   -1.0575 +102.5231 
  5   -1.2600 +101.7231 
  6   -1.4625 +100.9231 
  7   -1.6650 +100.1231 
  8   -1.8675  +99.3231 
  9   -2.0700  +98.5231 
 10   -2.2725  +97.7231 
 11   -2.4750  +96.9231 
 12   -2.6775  +96.1231 
 13   -2.8575  +95.3077 
 14   -3.0600  +94.5154 
 15   -3.2625  +93.7231 
 16   -3.4650  +92.9308 
 17   -3.6675  +92.1385 
 18   -3.8700  +91.3462 
 19   -4.0725  +90.5538 
 20   -4.2750  +89.7615 
 21   -4.4877  +88.9692 
 22   -4.6575  +88.1769 
 23   -4.8825  +87.3615 
 24   -5.0850  +86.5154 
 25   -5.2650  +85.7000 
 26   -5.4675  +84.9077 
 27   -5.6700  +84.1154 
 28   -5.8725  +83.3231 
 29   -6.0750  +82.5308 
 30   -6.2775  +81.7000 
 31   -6.5025  +80.8462 
 32   -6.6825  +80.0462 
 33   -6.8850  +79.2538 
 34   -7.0875  +78.4615 
 35   -7.2900  +77.6538 
 36   -7.5159  +76.7692 
 37   -7.6725  +75.9769 
 38   -7.8750  +75.1846 
 39   -8.1049  +74.3692 
 40   -8.2800  +73.5000 
 41   -8.4825  +72.7077 
 42   -8.6850  +71.9154 
 43   -8.9100  +71.0308 
 44   -9.0900  +70.2231 
 45   -9.2925  +69.4308 
 46   -9.5175  +68.5462 
 47   -9.6975  +67.7462 
 48   -9.9000  +66.9462 
 49  -10.1025  +66.0615 
 50  -10.3148  +65.2692 
 51  -10.4850  +64.3769 
 52  -10.6875  +63.5846 
 53  -10.9125  +62.7462 
 54  -11.0925  +61.9077 
 55  -11.2950  +61.0846 
 56  -11.4975  +60.2231 
 57  -11.7000  +59.3923 
 58  -11.9025  +58.5308 
 59  -12.1288  +57.6692 
 60  -12.3075  +56.8385 
 61  -12.5100  +55.9462 
 62  -12.7125  +55.1538 
 63  -12.9150  +54.2615 
 64  -13.1175  +53.4000 
 65  -13.2975  +52.5769 
 66  -13.5000  +51.6846 
 67  -13.7025  +50.7923 
 68  -13.9050  +50.0000 
 69  -14.1075  +49.1077 
 70  -14.3100  +48.2154 
 71  -14.5350  +47.3615 
 72  -14.7150  +46.5308 
 73  -14.9175  +45.6385 
 74  -15.1200  +44.7462 
 75  -15.3225  +43.8538 
 76  -15.5250  +42.9615 
 77  -15.7490  +42.0692 
 78  -15.9075  +41.2769 
 79  -16.1100  +40.3846 
 80  -16.3125  +39.4923 
 81  -16.5150  +38.6000 
 82  -16.7175  +37.7077 
 83  -16.9200  +36.8154 
 84  -17.1225  +35.9231 
 85  -17.3250  +34.9308 
 86  -17.5275  +34.0385 
 87  -17.7300  +33.1462 
 88  -17.9325  +32.2538 
 89  -18.1350  +31.3615 
 90  -18.3405  +30.4692 
 91  -18.5175  +29.4769 
 92  -18.7200  +28.5846 
 93  -18.9225  +27.6923 
 94  -19.1250  +26.8000 
 95  -19.3275  +25.8077 
 96  -19.5300  +24.9154 
 97  -19.7325  +23.9231 
 98  -19.9350  +23.0308 
 99  -20.1375  +22.1385 

[编辑3] 进展更新

对@Ben的一些澄清

工作原理

第一张图片下面的彩色方程式给出了半径r0，它由2个连接的90度三角形构成（基本三角法）。

红色内容：

• y(t)是电机位置，已知
• a(t)是IRC状态，也已知

绿色内容：

• a0,y0,z0是机械尺寸，已知但不精确，所以我测量了许多不同y(t)位置下的a(t)并使用已知校准管r0计算了更高精度的a0,y0,z0

进一步提高精度

我实际上通过使用特殊校准动作测量y1=y0+z0*cos(a0)来使其更加精确，精度约为0.03毫米。它是臂在a0位置和管道运动轴的交点高度，通过从上向下运动时臂第一次接触管道的情况来测量和插值出来的。它是由使用半径和a0重新计算的实际位置... 因为接触点不在这个轴上... (除非r0=0.0)。这也消除了校准中的一个近似循环，因为y1,a0,z0是相互依赖的，可以从彼此计算出来。另外，由于IRC的不连续测量方式以及a(t),y(t)位置的双重别名，消除了机器精度和计算稳定性方面的很多问题（在真正的机器上）。我现在无法可靠地评估精度，因为通过分析许多测量周期，我发现机器上存在一些机械问题，所以等待修复。无论如何，考虑到两种方法和_irc_calib_n=30，r0=80.03毫米的校准与模拟精度现在为：

    ;      computed     simulated  |delta|
    a0=  -6.915840 ;  -6.916710   +0.000870 deg
    y0=+186.009765 ;+186.012822   +0.003057 mm
    y1=+158.342452 ;+158.342187   +0.000264 mm
    z0=+228.102470 ;+228.100000   +0.002470 mm

较大的校准r0会降低精度（由于更有限的a(t)范围），这是通过计算所有a0,y0,(y1),z1得出的，没有直接测量或已知。虽然这已经可以接受，但如我之前所写，需要在机器准备就绪时进行检查。为了完整起见，现在模拟测量的样子如下图所示：

[编辑4]请参阅逼近搜索的工作原理

- Spektre

---

+1 非常详细的问题。不知道这是不是作业，但它在自己的领域里肯定很美。 - duffymo

@duffymo 不，这是我在工作中一直面临的一个问题....参数a0，y0，z0随时间变化，并且直接在机器上测量它们是不可行的，因此我寻找其他解决方案，而这是最接近我所需的。 - Spektre

@Ben，请查看我在问题中添加的一些澄清和进展更新… - Spektre

这似乎更适合在math.stackexchange.com上讨论。 - Alnitak

@Alnitak可能是这样，但数学专业人士通常不知道这些东西，我得到的最好的改进是通过改变运动模式并从不同的角度计算事物，这不是非常数学化的方法...方程式是正确的，所以没有太多可以改进的，由于导数的形状不好，迭代方法不可行...所以唯一的方法是最小误差搜索。但我仍然认为这可以更好地改进（不是从数学方面）。 - Spektre

Answer 1

如果我理解正确，您正在尝试通过测量y和a的值 而不是直接测量半径r0。

将常规误差传播应用于您对r0的公式，可以得到（对）r0的误差估计。在小角度极限下（适用于此处，因为a（t）限制为20度），这大致给出（使用三角函数的小角近似）

dr0^2 ~= dy^2 + z0^2 (pi*da/180)^2

因此，在r0远小于z0的情况下，r0的相对误差始终比y和z0*sin(a)的相对误差要大得多。这已经可以从您的图表中清楚地看出：测量量仅在一定程度上依赖于r0。

换句话说，这不是确定半径r0的聪明方法。除了增加角度a的范围之外，您无法解决这个根本限制。进行多次测量（以减少噪音/误差的常规方法）可能没有帮助，因为这些测量由于机器的内部工作而不是彼此独立的。因此，唯一的帮助就是更准确的测量。

为了分析情况，我建议制作图表/图形，例如推断的r0作为y的函数或y作为a的函数，以固定r0。

- Walter

---

你把它搞反了，r0 的计算是正确的（比我需要的更精确，半径的依赖关系已经足够了，但是这用于测量 r0），问题在于如何获得机械常数 a0、y0、z0 的精确值而不实际测量它们。因此，我在整个范围内测量曲线以了解半径（每个红色曲线在最后一张图像中都是单个半径的测量曲线），然后找到产生与测量曲线最接近的 a0、y0、z0，并将结果用于代替不准确测量的机械常数。在我需要的精度下无法进行测量... - Spektre

你从不会“测量”r0。把这个概念从你的脑海中排除掉。 - Walter

好的，对于 y0，您可以设置 r0=0，然后绘制 a(t) 相对于 y(t) 的图形并寻找极值点。至于精度，您可以考虑使用非常大的管子进行校准，其中 da/dy 变得非常大。 - Ben

在校准过程中，使用已知的测量半径r0，该半径可以被精确测量...在正常操作期间，半径是从a(t)和y(t)计算出来的...它已经实现并且正常工作，所以不要试图说服我（如果这不是真的，机器将无法正常运行）。问题在于机械常数a0，y0，z0可能会随时间改变，并且物理上测量它们是不可行的（这意味着必须拆卸机器并使用非常昂贵的测量设备，在实验室外不可能进行）。这个问题是关于校准而不是r0测量的... - Spektre

您IP地址为143.198.54.68，由于运营成本限制，当前对于免费用户的使用频率限制为每个IP每72小时10次对话，如需解除限制，请点击左下角设置图标按钮（手机用户先点击左上角菜单按钮）。 - Spektre