好吧,我玩得很开心。
创建一个名为 FuncEval 的类:
Option Explicit
Dim output_ As Double
Dim input_() As Double
Public Property Get VectArr() As Double()
VectArr = input_
End Property
Public Function Vect(i As Integer)
Vect = input_(i)
End Function
Public Sub SetVect(ByRef newVect() As Double)
Dim i As Integer
ReDim input_(LBound(newVect) To UBound(newVect)) As Double
For i = LBound(newVect) To UBound(newVect)
input_(i) = newVect(i)
Next i
End Sub
Public Property Get Result() As Double
Result = output_
End Property
Public Property Let Result(newRes As Double)
output_ = newRes
End Property
还有一个名为 Func 的类:
Option Explicit
Private cube_ As Double
Public Property Let Cube(newCube As Double)
cube_ = newCube
End Property
Public Function Eval(ByRef val() As Double) As FuncEval
Dim ret As New FuncEval
ret.Result = Abs(cube_ - val(0) * val(0) * val(0))
ret.SetVect val
Set Eval = ret
End Function
现在将此代码放入标准模块中:
Option Explicit
Function NelderMead(f As Func, _
ByRef guess() As Double) As Double()
'Algorithm follows that outlined here:
'http://www.mathworks.com/help/techdoc/math/bsotu2d.html#bsgpq6p-11
'Used as the perturbation for the initial guess when guess(i) == 0
Dim zeroPert As Double
zeroPert = 0.00025
'The factor each element of guess(i) is multiplied by to obtain the
'initial simplex
Dim pertFact As Double
pertFact = 1.05
'Tolerance
Dim eps As Double
eps = 0.000000000001
Dim shrink As Boolean
Dim i As Integer, j As Integer, n As Integer
Dim simplex() As Variant
Dim origVal As Double, lowest As Double
Dim m() As Double, r() As Double, s() As Double, c() As Double, cc() As Double, diff() As Double
Dim FE As FuncEval, FR As FuncEval, FS As FuncEval, FC As FuncEval, FCC As FuncEval, newFE As FuncEval
n = UBound(guess) - LBound(guess) + 1
ReDim m(LBound(guess) To UBound(guess)) As Double
ReDim r(LBound(guess) To UBound(guess)) As Double
ReDim s(LBound(guess) To UBound(guess)) As Double
ReDim c(LBound(guess) To UBound(guess)) As Double
ReDim cc(LBound(guess) To UBound(guess)) As Double
ReDim diff(LBound(guess) To UBound(guess)) As Double
ReDim simplex(LBound(guess) To UBound(guess) + 1) As Variant
Set simplex(LBound(simplex)) = f.Eval(guess)
'Generate the simplex
For i = LBound(guess) To UBound(guess)
origVal = guess(i)
If origVal = 0 Then
guess(i) = zeroPert
Else
guess(i) = pertFact * origVal
End If
Set simplex(LBound(simplex) + i - LBound(guess) + 1) = f.Eval(guess)
guess(i) = origVal
Next i
'Sort the simplex by f(x)
For i = LBound(simplex) To UBound(simplex) - 1
For j = i + 1 To UBound(simplex)
If simplex(i).Result > simplex(j).Result Then
Set FE = simplex(i)
Set simplex(i) = simplex(j)
Set simplex(j) = FE
End If
Next j
Next i
Do
Set newFE = Nothing
shrink = False
lowest = simplex(LBound(simplex)).Result
'Calculate m
For i = LBound(m) To UBound(m)
m(i) = 0
For j = LBound(simplex) To UBound(simplex) - 1
m(i) = m(i) + simplex(j).Vect(i)
Next j
m(i) = m(i) / n
Next i
'Calculate the reflected point
For i = LBound(r) To UBound(r)
r(i) = 2 * m(i) - simplex(UBound(simplex)).Vect(i)
Next i
Set FR = f.Eval(r)
'Check acceptance conditions
If (simplex(LBound(simplex)).Result <= FR.Result) And (FR.Result < simplex(UBound(simplex) - 1).Result) Then
'Accept r, replace the worst value and iterate
Set newFE = FR
ElseIf FR.Result < simplex(LBound(simplex)).Result Then
'Calculate the expansion point, s
For i = LBound(s) To UBound(s)
s(i) = m(i) + 2 * (m(i) - simplex(UBound(simplex)).Vect(i))
Next i
Set FS = f.Eval(s)
If FS.Result < FR.Result Then
Set newFE = FS
Else
Set newFE = FR
End If
ElseIf FR.Result >= simplex(UBound(simplex) - 1).Result Then
'Perform a contraction between m and the better of x(n+1) and r
If FR.Result < simplex(UBound(simplex)).Result Then
'Contract outside
For i = LBound(c) To UBound(c)
c(i) = m(i) + (r(i) - m(i)) / 2
Next i
Set FC = f.Eval(c)
If FC.Result < FR.Result Then
Set newFE = FC
Else
shrink = True
End If
Else
'Contract inside
For i = LBound(cc) To UBound(cc)
cc(i) = m(i) + (simplex(UBound(simplex)).Vect(i) - m(i)) / 2
Next i
Set FCC = f.Eval(cc)
If FCC.Result < simplex(UBound(simplex)).Result Then
Set newFE = FCC
Else
shrink = True
End If
End If
End If
'Shrink if required
If shrink Then
For i = LBound(simplex) + 1 To UBound(simplex)
For j = LBound(simplex(i).VectArr) To UBound(simplex(i).VectArr)
diff(j) = simplex(LBound(simplex)).Vect(j) + (simplex(i).Vect(j) - simplex(LBound(simplex)).Vect(j)) / 2
Next j
Set simplex(i) = f.Eval(diff)
Next i
End If
'Insert the new element in place
If Not newFE Is Nothing Then
For i = LBound(simplex) To UBound(simplex)
If simplex(i).Result > newFE.Result Then
For j = UBound(simplex) To i + 1 Step -1
Set simplex(j) = simplex(j - 1)
Next j
Set simplex(i) = newFE
Exit For
End If
Next i
End If
Loop Until (simplex(UBound(simplex)).Result - simplex(LBound(simplex)).Result) < eps
NelderMead = simplex(LBound(simplex)).VectArr
End Function
Function test(cube, guess) As Double
Dim f As New Func
Dim guessVec(0 To 0) As Double
Dim Result() As Double
Dim i As Integer
Dim output As String
f.cube = cube
guessVec(0) = guess
Result = NelderMead(f, guessVec)
test = Result(0)
End Function
Func 类包含您的残差函数。NelderMead 方法只需要 Func 类的 Result 方法,因此只要 Eval 方法处理与您的初始猜测长度相同的向量并返回 FuncEval 对象,您就可以对 Func 类做任何事情。
调用测试函数来查看它的运行情况。注意,我还没有实际测试过多维向量,我必须出去,如果您有任何问题,请告诉我!
编辑:关于泛化函数传递的建议
您需要针对不同的问题创建许多不同的类。这意味着要保持 NelderMead 函数的通用性,您需要将其声明行更改为以下内容:
Function NelderMead(f As Object, _
ByRef guess() As Double) As Double()
无论 f 是什么,它都必须始终有一个 Eval 方法,该方法采用一个双精度数组。
编辑:函数传递,可能是在 VBA 中完成的(愚蠢的)方式
Function f(x() As Double) As Double
f = x(0) * x(0)
End Function
Sub Test()
Dim x(0 To 0) As Double
x(0) = 5
Debug.Print Application.Run("f", x)
End Sub
使用此方法,您将拥有以下声明:
Function NelderMead(f As String, _
ByRef guess() As Double) As Double()
然后使用上面的 Application.Run 语法调用 f 。您还需要在函数内部进行一些更改。它并不漂亮,但坦率地说,它一开始并不那么漂亮。