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
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
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
下面用求组合各个圆的斜率的平均值作为最终的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
发现这样求平均和线性规划差别挺大的。所以还是用线性规划吧。
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
上面代码比较适合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