最小二乘法的应用-C语言实现线损电压补偿-代码篇3

最小二乘法的应用-C语言实现线损电压补偿-求出最小二乘法的系数

最小二乘法的基本原理

原理篇1

一段不可精确测出的回路电阻

应用场景2

求出最小二乘法的系数

表达式

1. 数学表达式

2. 转化过程
   最小二乘法求回归直线方程的推导过程

式1
sumX * sumY - n * sumXY = (n * meanX) * (n * meanY) - n * sigama^(Xi * Yi) =
-n * (sigama^(Xi * Yi) - n * meanX * meanY)

式2
sumX * sumX - n * sumXX = (n * meanX) * (n * meanX) - n * sumXX =
-n * (sigama^(Xi * Xi) - n * meanX * meanX)

b = 式1 / 式2 =
(sumX * sumY - n * sumXY) / (sumX * sumX - n * sumXX)

式3
(sumY - b * sumX) / n = (n * meanY - b * n * meanX) / n =
meanY - b * meanX

a = 式3

3. MCU单片机应用表达式
   使用移位的方法, 使精度保留到一定的位数, 减少失真;

y = ( norm_slope * x + norm_offset ) >> norm
norm_slope = ( ( sumX * sumY - n * sumXY ) << norm ) / ( sumX * sumX - n * sumXX )
norm_offset = ( ( sumY << norm ) - norm_slope * sumX ) / n

OK. Talk is cheap. Show you the code.

斜率计算方法(贪方便用C#去读写文件处理数据)

        private void LeastSquareMethod(in int[] intArrX, in int[] intArrY, in int numOfArr, out int sdwSlope, out int sdwOffset, out int wNorm)
        {
            Int32 sdwTemp = 0;
            UInt32 sumX = 0, sumY = 0, sumXY = 0, sumXX = 0;
            for (int i = 0; i < numOfArr; i++)
            {
                UInt32 x = 0, y = 0;
                x = (UInt32)intArrX[i];
                sumX += x;
                sumXX += x * x;

                y = (UInt32)intArrY[i];
                sumY += y;
                sumXY += x * y;
            }

            try
            {
                checked
                {
                    // Calculate the normalized slope and offset using the following equations
                    //      y = ( norm_slope * x + norm_offset ) >> norm
                    //      norm_slope = ( ( sumX * sumY - n * sumXY ) << norm ) / ( sumX * sumX - n * sumXX )
                    //      norm_offset = ( ( sumY << norm ) - norm_slope * sumX ) / n

                    // Calculate the denominator of the coil current slope
                    sdwTemp = ((Int32)sumX * (Int32)sumX - (Int32)numOfArr * (Int32)sumXX);

                    // Calculate the numerator of the coil current slope 
                    sdwSlope = ((Int32)sumX * (Int32)sumY - (Int32)numOfArr * (Int32)sumXY);

                    // Increment the normalization value until the slope is about to overrun
                    // Note: 1073741824 = 0x40000000
                    wNorm = 0;
                    while ((sdwSlope != 0) && (sdwSlope < 1073741824) && (sdwSlope > -1073741824))
                    {
                        sdwSlope = (sdwSlope * 2);
                        wNorm++;
                    }

                    sdwSlope /= sdwTemp;

                    // Calculate the offset
                    sdwOffset = (((Int32)sumY * (Int32)(1UL << wNorm)) - sdwSlope * (Int32)sumX) / (Int32)numOfArr;

                    // Save the slope and offset values to the calibration structure
                    //pCoilCurFrmRailVol->sdwSlope = sdwSlope;
                    //pCoilCurFrmRailVol->sdwOffset = sdwOffset;
                    //pCoilCurFrmRailVol->wNorm = wNorm;

                    Console.WriteLine("sdwSlope = {0}, sdwOffset = {1}, wNorm = {2}", sdwSlope, sdwOffset, wNorm);

                }
            }
            catch (OverflowException ex)
            {
                sdwSlope = 0;
                sdwOffset = 0;
                wNorm = 0;
                MessageBox.Show(ex.Message + "\n亦即 OverFlow，" +
                    "\n请检查选中目标文件数据是否正确，" +
                    "\n以及数据的数值是否太大，可否考虑略去数值较大的数据?");
            }
            catch (Exception e)
            {
                sdwSlope = 0;
                sdwOffset = 0;
                wNorm = 0;
                MessageBox.Show(e.Message);
            }
        }

最终得到的系数为sdwSlope(斜率b) sdwOffset(偏移值a) wNorm(移位保精度)

MCU计算时调用方法

这里, 回路线损deltaVbat = (wInputPower * sdwSlope + sdwOffset) >> wNorm;

#if VBAT_COMPENSATE_BASED_ON_INPUT_POWER
uint16 append_PROT_CheckBoardParams_GetDeltaVbatJudgeByInputPower (uint16 wInputPower)
{
	sint32 sdwTemp;
	uint32  dwTemp;
    CURRENT_FITTING_PARAMS DeltaVbatFrmCurrent;
	CURRENT_FITTING_PARAMS *pDVbatFrmCur = &DeltaVbatFrmCurrent;
    pDVbatFrmCur->sdwSlope  = VBAT_COMPENSATE_BASED_ON_INPUT_POWER_SLOPE;
    pDVbatFrmCur->sdwOffset = VBAT_COMPENSATE_BASED_ON_INPUT_POWER_OFFSET;
    pDVbatFrmCur->wNorm     = VBAT_COMPENSATE_BASED_ON_INPUT_POWER_NORM;
    
	sdwTemp = wInputPower / 10;
	sdwTemp *= pDVbatFrmCur->sdwSlope;
	sdwTemp += pDVbatFrmCur->sdwOffset;
	if(sdwTemp < 0)
	{
	    //in case of error
		sdwTemp = 0;
	}
	dwTemp = sdwTemp;
	dwTemp >>= pDVbatFrmCur->wNorm;
	return (uint16)dwTemp;
}
#endif // EOF #if VBAT_COMPENSATE_BASED_ON_INPUT_POWER

好啦. 就甘啦. 计算结果呢, 在误差要求范围内, 甲方爸爸 客户 满意就行;