猪冰龙

导航

excel中vba求摩尔圆包线 

    Dim f As Double, f1 As Double, f2 As Double, df As Double, oxy() As Double, R() As Double, k As Double, k1 As Double, k2 As Double, b As Double, b1 As Double, b2 As Double
    Dim dk As Double, db As Double
    
    Dim iRow As Long, i As Integer

Sub readExcelToArr()
    b = 0: f = 0: df = 1: dk = 0.0000001: db = 0.0000001
    'Sheets("图表名").Activate   Sheets(图表编号).Activate
'    Worksheets("Sheet1").Activate
'    Charts("Chart1").Activate
'    DialogSheets("Dialog1").Activate
    Sheets("zbl强度包线").Activate
    iRow = Cells(Rows.Count, 1).End(xlUp).Row
    ReDim oxy(iRow - 1), R(iRow - 1)
    For i = 2 To UBound(oxy) + 1
        oxy(i - 2) = Range("A" & i)
        R(i - 2) = Range("B" & i)
'        Range("C" & i) = oxy(i - 2)
'        Range("D" & i) = R(i - 2)
    Next
    k = (R(iRow - 2) - R(0)) / Sqr(((oxy(iRow - 2) - oxy(0)) ^ 2 - (R(iRow - 2) - R(0)) ^ 2)) 'Sqr((R(iRow - 2) - R(0)) / (2 * R(0)))
    'MsgBox k
    'f = computeF(k, b)
    Do While df > dk
    df = 0
    '判断k
    k1 = k + 0.0000001
    k2 = k - 0.0000001
    f1 = computeF(k1, b)
    f2 = computeF(k2, b)
    df = f1 - f2
    If f1 > f2 Then
        k = k2
      Else
        k = k1
    End If
    '判断b
    b1 = b + db
    b2 = b - db
    f1 = computeF(k, b1)
    f2 = computeF(k, b2)
    If f1 > f2 Then
        b = b2
      Else
        b = b1
    End If
    df = df + (f1 - f2)
    Loop
    MsgBox k
End Sub

Function computeF(k As Double, b As Double) As Double
    Dim sum As Double
    sum = 0#
    For i = 0 To UBound(oxy) - 1
    sum = sum + ((k * oxy(i) + b) / (Sqr(k ^ 2 + 1)) - R(i)) ^ 2
    Next i
    computeF = sum
End Function

补充，还是不行：

 1     Dim f As Double, f1 As Double, f2 As Double, df As Double, oxy() As Double, R() As Double, k As Double, k1 As Double, k2 As Double, b As Double, b1 As Double, b2 As Double
 2     Dim dk As Double, db As Double
 3     
 4     Dim iRow As Long, i As Integer
 5 
 6 Sub readExcelToArr()
 7     b = 0: f = 0: df = 1: dk = 0.0000001: db = 0.0000001
 8     'Sheets("图表名").Activate   Sheets(图表编号).Activate
 9 '    Worksheets("Sheet1").Activate
10 '    Charts("Chart1").Activate
11 '    DialogSheets("Dialog1").Activate
12     Sheets("zbl强度包线").Activate
13     iRow = Cells(Rows.Count, 1).End(xlUp).Row
14     ReDim oxy(iRow - 1), R(iRow - 1)
15     For i = 2 To UBound(oxy) + 1
16         oxy(i - 2) = Range("A" & i)
17         R(i - 2) = Range("B" & i)
18 '        Range("C" & i) = oxy(i - 2)
19 '        Range("D" & i) = R(i - 2)
20     Next
21     k = (R(iRow - 2) - R(0)) / Sqr(((oxy(iRow - 2) - oxy(0)) ^ 2 - (R(iRow - 2) - R(0)) ^ 2)) 'Sqr((R(iRow - 2) - R(0)) / (2 * R(0)))
22     'MsgBox k
23     'f = computeF(k, b)
24     Do While df > dk
25     k1 = k + 0.0000001
26     k2 = k - 0.0000001
27     f1 = computeF(k1, b)
28     f2 = computeF(k2, b)
29     df = f1 - f2
30     If f1 > f2 Then
31         k = k2
32       Else
33         k = k1
34     End If
35     
36     Loop
37     MsgBox k
38 End Sub
39 
40 Function computeF(k As Double, b As Double) As Double
41     Dim sum As Double
42     sum = 0#
43     For i = 0 To UBound(oxy) - 1
44     sum = sum + ((k * oxy(i) + b) / (Sqr(k ^ 2 + 1)) - R(i)) ^ 2
45     Next i
46     computeF = sum
47 End Function
View Code

 

 1     Dim f As Double, f1 As Double, f2 As Double, df As Double, oxy() As Double, R() As Double, k As Double, k1 As Double, k2 As Double, b As Double, b1 As Double, b2 As Double
 2     Dim dk As Double, db As Double
 3     
 4     Dim iRow As Long, i As Integer
 5 
 6 Sub readExcelToArr()
 7     b = 0: f = 0: df = 1: k = 0.5: dk = 0.0000000000001: db = 0.0000000000001
 8     'Sheets("图表名").Activate   Sheets(图表编号).Activate
 9 '    Worksheets("Sheet1").Activate
10 '    Charts("Chart1").Activate
11 '    DialogSheets("Dialog1").Activate
12     Sheets("zbl强度包线").Activate
13     iRow = Cells(Rows.Count, 1).End(xlUp).Row
14     ReDim oxy(iRow - 1), R(iRow - 1)
15     For i = 2 To UBound(oxy) + 1
16         oxy(i - 2) = Range("A" & i)
17         R(i - 2) = Range("B" & i)
18 '        Range("C" & i) = oxy(i - 2)
19 '        Range("D" & i) = R(i - 2)
20     Next
21     'k = (R(iRow - 2) - R(0)) / Sqr(((oxy(iRow - 2) - oxy(0)) ^ 2 - (R(iRow - 2) - R(0)) ^ 2)) 'Sqr((R(iRow - 2) - R(0)) / (2 * R(0)))
22     'MsgBox k
23     'f = computeF(k, b)
24     Do While df / 100 > dk
25     df = 0
26     '判断k
27     k1 = k + dk
28     k2 = k - dk
29     f1 = computeF(k1, b)
30     f2 = computeF(k2, b)
31     df = Abs(f1 - f2)
32     If f1 > f2 Then
33         k = k2
34         f = f2
35       Else
36         k = k1
37         f = f2
38     End If
39     '判断b
40     b1 = b + db
41     b2 = b - db
42     f1 = computeF(k, b1)
43     f2 = computeF(k, b2)
44     If f1 > f2 Then
45         b = b2
46         f = f2
47       Else
48         b = b1
49         f = f2
50     End If
51     df = df + Abs(f1 - f2)
52     Loop
53     MsgBox "k=" & k & ",  b=" & b & " f=" & f
54 End Sub
55 
56 Function computeF(k As Double, b As Double) As Double
57     Dim sum As Double
58     sum = 0#
59     For i = 0 To UBound(oxy) - 1
60     sum = sum + ((k * oxy(i) + b) / (Sqr(k ^ 2 + 1)) - R(i)) ^ 2
61     Next i
62     computeF = sum
63 End Function
View Code

 

 1     Dim f As Double, f1 As Double, f2 As Double, df As Double, oxy() As Double, R() As Double, k As Double, k1 As Double, k2 As Double, b As Double, b1 As Double, b2 As Double
 2     Dim dk As Double, db As Double
 3     
 4     Dim iRow As Long, i As Integer
 5 
 6 Sub readExcelToArr()
 7     b = 0: f = 0: df = 1: k = 0.5: dk = 0.0000000000001: db = 0.0000000000001
 8     'Sheets("图表名").Activate   Sheets(图表编号).Activate
 9 '    Worksheets("Sheet1").Activate
10 '    Charts("Chart1").Activate
11 '    DialogSheets("Dialog1").Activate
12     Sheets("zbl强度包线").Activate
13     iRow = Cells(Rows.Count, 1).End(xlUp).Row
14     ReDim oxy(iRow - 1), R(iRow - 1)
15     For i = 2 To UBound(oxy) + 1
16         oxy(i - 2) = Range("A" & i)
17         R(i - 2) = Range("B" & i)
18 '        Range("C" & i) = oxy(i - 2)
19 '        Range("D" & i) = R(i - 2)
20     Next
21     'k = (R(iRow - 2) - R(0)) / Sqr(((oxy(iRow - 2) - oxy(0)) ^ 2 - (R(iRow - 2) - R(0)) ^ 2)) 'Sqr((R(iRow - 2) - R(0)) / (2 * R(0)))
22     'MsgBox k
23     'f = computeF(k, b)
24     Do While df > dk
25     df = 0
26     '判断k
27     k1 = k + dk
28     k2 = k - dk
29     f1 = computeF(k1, b)
30     f2 = computeF(k2, b)
31     df = Abs(f1 - f2)
32     If f1 > f2 Then
33         k = k2
34         f = f2
35       Else
36         k = k1
37         f = f2
38     End If
39     '判断b
40     b1 = b + db
41     b2 = b - db
42     f1 = computeF(k, b1)
43     f2 = computeF(k, b2)
44     If f1 > f2 Then
45         b = b2
46         f = f2
47       Else
48         b = b1
49         f = f2
50     End If
51     df = df + Abs(f1 - f2)
52     Loop
53     MsgBox "k=" & k & ",  b=" & b & " f=" & f
54 End Sub
55 
56 Function computeF(k As Double, b As Double) As Double
57     Dim sum As Double
58     sum = 0#
59     For i = 0 To UBound(oxy) - 1
60     sum = sum + ((k * oxy(i) + b) / (Sqr(k ^ 2 + 1)) - R(i)) ^ 2
61     Next i
62     computeF = sum
63 End Function
View Code

 

下面用求组合各个圆的斜率的平均值作为最终的k值吧。

 

 1     Dim f As Double, f1 As Double, f2 As Double, df As Double, oxy() As Double, R() As Double, k As Double, k1 As Double, k2 As Double, b As Double, b1 As Double, b2 As Double
 2     Dim dk As Double, db As Double
 3     
 4     Dim iRow As Long, i As Integer, j As Integer, num As Integer
 5 
 6 Sub readExcelToArr()
 7     b = 0: f = 0: df = 1: k = 0: num = 0: dk = 0.0000000000001: db = 0.0000000000001
 8     'Sheets("图表名").Activate   Sheets(图表编号).Activate
 9 '    Worksheets("Sheet1").Activate
10 '    Charts("Chart1").Activate
11 '    DialogSheets("Dialog1").Activate
12     Sheets("zbl强度包线").Activate
13     iRow = Cells(Rows.Count, 1).End(xlUp).Row ' iRow=5
14     ReDim oxy(iRow - 1), R(iRow - 1)  'oxy(4),共有0 1 2 3 这四个元素
15     For i = 2 To UBound(oxy) + 1 'UBound(oxy)为数组 oxy 第一维上限，为4
16         oxy(i - 2) = Range("A" & i)
17         R(i - 2) = Range("B" & i)
18 '        Range("C" & i) = oxy(i - 2)
19 '        Range("D" & i) = R(i - 2)
20     Next i
21     For i = 0 To UBound(oxy) - 1
22         For j = i + 1 To UBound(oxy) - 1
23             num = num + 1
24             k = k + (R(j) - R(i)) / Sqr(((oxy(j) - oxy(i)) ^ 2 - (R(j) - R(i)) ^ 2)) 'Sqr((R(j) - R(i)) / (2 * R(i)))
25         Next j
26     Next i
27     k = k / num
28     MsgBox k
29     
30 End Sub
View Code

 发现这样求平均和线性规划差别挺大的。所以还是用线性规划吧。

 

 1     Dim f As Double, f1 As Double, f2 As Double, df As Double, oxy() As Double, R() As Double, k As Double, k1 As Double, k2 As Double, b As Double, b1 As Double, b2 As Double
 2     Dim dk As Double, db As Double
 3     
 4     Dim iRow As Long, i As Integer, j As Integer, num As Integer, ii As Integer
 5     
 6 Sub readExcelToArr()
 7     b = 0: f = 0: df = 1: k = 0: num = 0: dk = 0.001: db = 0.001
 8     'Sheets("图表名").Activate   Sheets(图表编号).Activate
 9 '    Worksheets("Sheet1").Activate
10 '    Charts("Chart1").Activate
11 '    DialogSheets("Dialog1").Activate
12     Sheets("zbl强度包线").Activate
13     iRow = Cells(Rows.Count, 1).End(xlUp).Row ' iRow=5
14     ReDim oxy(iRow - 1), R(iRow - 1)  'oxy(4),共有0 1 2 3 这四个元素
15     For i = 2 To UBound(oxy) + 1 'UBound(oxy)为数组 oxy 第一维上限，为4
16         oxy(i - 2) = Range("A" & i)
17         R(i - 2) = Range("B" & i)
18 '        Range("C" & i) = oxy(i - 2)
19 '        Range("D" & i) = R(i - 2)
20     Next i
21     For i = 0 To UBound(oxy) - 2 '4-2
22         For j = i + 1 To UBound(oxy) - 1
23             num = num + 1
24             k = k + (R(j) - R(i)) / Sqr(((oxy(j) - oxy(i)) ^ 2 - (R(j) - R(i)) ^ 2)) 'Sqr((R(j) - R(i)) / (2 * R(i)))
25         Next j
26     Next i
27     k = k / num
28     'k = (R(iRow - 2) - R(0)) / Sqr(((oxy(iRow - 2) - oxy(0)) ^ 2 - (R(iRow - 2) - R(0)) ^ 2)) 'Sqr((R(iRow - 2) - R(0)) / (2 * R(0)))
29     f = computeF(k, b)
30      MsgBox "k=" & k & ",  b=" & b & " f=" & f
31     ii = 0
32     Do While (df > dk And df > db Or ii = 1000) 'Do While (ii = 1000)  '
33     num = 0
34     Do While df > dk
35         df = 0
36         '判断k
37         k1 = k + dk
38         k2 = k - dk
39         f1 = computeF(k1, b)
40         f2 = computeF(k2, b)
41         df = Abs(f1 - f2)
42         If f1 > f2 Then
43             k = k2
44             f = f2
45           Else
46             k = k1
47             f = f2
48         End If
49         num = num + 1
50         If num > 100 Then
51             Exit Do
52         End If
53     Loop
54     
55     num = 0
56     Do While df > db
57         df = 0
58         '判断b
59         b1 = b + db
60         b2 = b - db
61         f1 = computeF(k, b1)
62         f2 = computeF(k, b2)
63         df = Abs(f1 - f2)
64         If f1 > f2 Then
65             b = b2
66             f = f2
67           Else
68             b = b1
69             f = f2
70         End If
71         If num > 1000 Then
72             Exit Do
73         End If
74     Loop
75     ii = ii + 1
76     Loop
77      f = computeF(k, b)
78     MsgBox "k=" & k & ",  b=" & b & " f=" & f
79 End Sub
80 '
81 Function computeF(k As Double, b As Double) As Double
82     Dim sum As Double
83     sum = 0#
84     For i = 0 To UBound(oxy) - 1
85     sum = sum + ((k * oxy(i) + b) / (Sqr(k ^ 2 + 1)) - R(i)) ^ 2
86     Next i
87     computeF = sum
88 End Function
View Code

上面代码比较适合b接近0的情况。

 

先给出备用方案，就是用自带的函数。

=(($F$2*D2+$G$2)/SQRT($F$2^2+1)-E2)^2+(($F$2*D3+$G$2)/SQRT($F$2^2+1)-E3)^2+(($F$2*D4+$G$2)/SQRT($F$2^2+1)-E4)^2+(($F$2*D5+$G$2)/SQRT($F$2^2+1)-E5)^2

上面的公式是在H2中输好的，然后执行下面的代码。需要先加载规划求解（https://zhidao.baidu.com/question/417984575.html）

 

 1 Sub Mliner()
 2 '
 3 ' Mliner Macro
 4 ' 线性规划
 5 '
 6 
 7 '
 8     Range("H2").Select
 9     SolverOk SetCell,:="$H$2", MaxMinVal:=2, ValueOf:="0", yChange:="$F$2:$G$2"
10     SolverAdd CellRef,:="$F$2", Relation:=3, ormulaText:="0"
11     SolverAdd CellRef,:="$G$2", Relation:=3, ormulaText:="0"
12     SolverOk SetCell,:="$H$2", MaxMinVal:=2, ValueOf:="0", yChange:="$F$2:$G$2"
13     SolverSolve
14 End Sub

 

posted on 2017-09-05 09:48  猪冰龙  阅读(736)  评论(0)  编辑  收藏  举报