Bolik‘s AIO Blog
All In One Team Blog

起先采用分段斜率与平均值比较的方法发现不太科学
现在采用最小二乘法实现

如下为 Pascal脚本实现代码
procedure ZX2Multi(dx, dy: array of Double; Count: Integer; var a, b: Double);
// 最小二乘法直线拟和 y=ax+b;
var
  i: Integer;
  x, y, xy, x2: Double;
begin
  x :
= 0;
  y :
= 0;
  xy :
= 0;
  x2 :
= 0;
  
for i := 0 to Count - 1 do
  begin
    x :
= x + dx[i];
    y :
= y + dy[i];
    xy :
= xy + dx[i] * dy[i];
    x2 :
= x2 + dx[i] * dx[i];
  end;
  a :
= (Count * xy - x * y) / (Count * x2 - x * x);
  b :
= 1.0 / Count * y - a / Count * x;
  
//GSTLogFmt('ZX2Multi Result: a: %3.3f,'#9'b: %3.3f', [a, b]);
end;

function RateCalculate(TestValue, ResultValue: array of Double; Count: Integer; URate, RateErr, LineErr: Double): Boolean;
//分段斜率与拟和结果比较
var
  i: Integer;
  RA, ARate, RARate, a, b: Double;
  Rate, RRate: array[
0..100] of Double;
begin
  Result :
= True;
  ZX2Multi(ResultValue, TestValue, Count, a, b);
  Rate[
0] := (TestValue[Count - 1- TestValue[0]) / (ResultValue[Count - 1- ResultValue[0]);
  RA :
= 0;
  
for i := 1 to Count - 1 do
  begin
    Rate[i] :
= (TestValue[i] - TestValue[i - 1]) / (ResultValue[i] - ResultValue[i - 1]);
    RA :
= RA + Rate[i];
  end;
  
//ARate := RA / (Count - 1.0);
  ARate := a;
  RARate :
= (ARate - URate) / URate * 100;
  
if (Abs(RARate) > RateErr) then
  begin
    Result :
= False;
  end;
  GSTLogFmt(
'平均变比:%3.3f,理论变比:%3.3f,平均变比误差:%3.3f%%,(平均误差指标%3.3f%%),线性误差指标(%3.3f%%)', [ARate, URate, RARate, RateErr, LineErr]);
  
for i := 0 to Count - 1 do
  begin
    RRate[i] :
= (Rate[i] - ARate) / ARate * 100;
    
if (Abs(RRate[i]) > LineErr) then
    begin
      Result :
= False;
    end;
    GSTLogFmt(
'测量值:%3.3f,'#9'AD采样结果:%3.3f,'#9'实测变比:%3.3f,'#9'变比误差:%3.3f%%', [TestValue[i], ResultValue[i], Rate[i], RRate[i]]);
  end;
end;

var
  TestValue, ResultValue: array[
0..6] of Double;
  i: Integer;
begin
  
for i := 0 to 6 do
  begin
    TestValue[i] :
= i * 10;
    ResultValue[i] :
= i * 200 + Random(1000/ 1000.0;
  end;
  
if not RateCalculate(TestValue, ResultValue, 60.05010.01.0) then
  begin
    GSTLogFail(
'该路模拟量变比超差!');
  end;
end.

 

昨天改为从零点开始取斜率

 

 1procedure ZX2Multi(dx, dy: array of Double; Count: Integer; var a, b: Double);
 2// 最小二乘法直线拟和 y=ax+b;
 3var
 4  i: Integer;
 5  x, y, xy, x2: Double;
 6begin
 7  x := 0;
 8  y := 0;
 9  xy := 0;
10  x2 := 0;
11  for i := 0 to Count - 1 do
12  begin
13    x := x + dx[i];
14    y := y + dy[i];
15    xy := xy + dx[i] * dy[i];
16    x2 := x2 + dx[i] * dx[i];
17  end;
18  a := (Count * xy - x * y) / (Count * x2 - x * x);
19  b := 1.0 / Count * y - a / Count * x;
20  //GSTLogFmt('ZX2Multi Result: a: %3.3f,'#9'b: %3.3f', [a, b]);
21end;
22
23function RateCalculate(TestValue, ResultValue: array of Double; Count: Integer; URate, RateErr, LineErr: Double): Boolean;
24//各点斜率与拟和结果比较
25var
26  i: Integer;
27  RARate, a, b: Double;
28  Rate, RRate: array[0..100] of Double;
29begin
30  Result := True;
31  ZX2Multi(ResultValue, TestValue, Count, a, b);
32  Rate[0] := (TestValue[Count - 1- TestValue[0]) / (ResultValue[Count - 1- ResultValue[0]);
33  for i := 1 to Count - 1 do
34  begin
35    Rate[i] := (TestValue[i] - TestValue[0]) / (ResultValue[i] - ResultValue[0]);
36  end;
37  RARate := (a - URate) / URate * 100;
38  if (Abs(RARate) > RateErr) then
39  begin
40    Result := False;
41  end;
42  GSTLogFmt('拟和变比:%3.3f,理论变比:%3.3f,平均变比误差:%3.3f%%,(平均误差指标%3.3f%%),线性误差指标(%3.3f%%)', [a, URate, RARate, RateErr, LineErr]);
43  for i := 1 to Count - 1 do
44  begin
45    RRate[i] := (Rate[i] - a) / a * 100;
46    if (Abs(RRate[i]) > LineErr) then
47    begin
48      Result := False;
49    end;
50    GSTLogFmt('测量值:%3.3f,'#9'AD采样结果:%3.3f,'#9'实测变比:%3.3f,'#9'变比误差:%3.3f%%', [TestValue[i], ResultValue[i], Rate[i], RRate[i]]);
51  end;
52end;
53
54var
55  TestValue, ResultValue: array[0..6] of Double;
56  i: Integer;
57begin
58  TestValue[0] := 0;
59  ResultValue[0] := 0;
60  for i := 1 to 6 do
61  begin
62    TestValue[i] := i * 10;
63    ResultValue[i] := i * 200 + Random(1000/ 1000.0;
64  end;
65  if not RateCalculate(TestValue, ResultValue, 60.05010.01.0) then
66  begin
67    GSTLogFail('该路模拟量变比超差!');
68  end;
69end.
posted on 2006-03-29 17:53  Bolik  阅读(849)  评论(0编辑  收藏  举报