改良单纯形法的VBA实现

付永超 / 2019-03-30

今年年初，我是使用VBA实现改良的单纯形法，该算法是参考《最优化计算方法》，“第八章无约束优化方法，8.3可变单纯形法的实现”章节实现的，该书的作者是蒋金山。这本书步骤非常详细的说明了计算方法和步骤，非常适合学习与练习使用。

实现的目的当然不是为了重新造一个轮子，主要是为了让自己对该算法有进一步的认识。通过自己动手编程的练习，让自己了解了“单纯形”是什么怎么进行反射，怎么算质心，怎么迭代？等等。由此让自己真正的了解了该算法。

代码如下：

Option Explicit

Dim xa, xb, xr, xe, zongi '反射系数，压缩系数，扩展系数，允许误差，迭代总次数

Dim fg As Range, fl As Range, fh As Range '三组参数中的最大、中间、最小

Dim ps1 As Range, ps2 As Range, ps3 As Range '三组参数中的第一列、第二列、第三列
Dim nx As Range, n3x As Range, n4x As Range '储存中间结果的下一步，下下一步的参数值
Sub rukou()
'入口程序，声明一些关键的变量，并串联起来求解过程
'该算法是参考《最优化计算方法》，蒋金山，第八章无约束优化方法，8.3可变单纯形法的实现
'可变单纯形法是一种直接搜索方法，1962年被提出，1964年给改进
Dim pp
Dim a, b, c
Dim rowi, ls1
zongi = 0 '迭代总次数
xa = 1 '反射系数
xb = 0.5 '压缩系数
xr = 2 '扩展系数
xe = 0.0000000002 '允许误差

'清除上一次的迭代结果
ls1 = ActiveSheet.Cells(9999, 16).End(xlUp).Row + 1
ActiveSheet.Range(ActiveSheet.Cells(2, 15), ActiveSheet.Cells(ls1, 20)).Clear


Do
    Call duquhuanjing2 '读取表格中给如入的3组初始的参数值，存入ps1，ps2，ps3变量中去
    Call d1_px(ps1, ps2, ps3) '为3组初始参数进行排序fh,fg,gl，以便接下来的计算，直接使用其中的最大值与最小值
    Call d2_mean(fh, fg, fl) '应用3组参数，计算新的更好的点替代fh
    pp = panju(fh, fg, fl) '第8步，收敛标准值的计算
    
    zongi = zongi + 1 '迭代总次数
    a = fh(1) '当前fh的目标函数值
    b = fg(1) '当前fg的目标函数值
    c = fl(1) '当前fl的目标函数值
    rowi = ActiveSheet.Cells(9999, 16).End(xlUp).Row + 1 '获取迭代结果显示的最后一行之后的一行行号
    ActiveSheet.Cells(rowi, 20) = pp '收敛标准的值
    Debug.Print zongi; "::fh, fg, fl="; a; b; c
    '将迭代结果输出到Excel表格中
    ActiveSheet.Cells(rowi, 16) = zongi
    ActiveSheet.Cells(rowi, 17) = fl(0)
    ActiveSheet.Cells(rowi, 18) = fg(0)
    ActiveSheet.Cells(rowi, 19) = fh(0)
   
    '清空操作对象的内容，便于重新开始迭代
    Set fh = Nothing
    Set fg = Nothing
    Set fl = Nothing
    Set ps1 = Nothing
    Set ps2 = Nothing
    Set ps3 = Nothing
    Set nx = Nothing
    Set n3x = Nothing
    Set n4x = Nothing
    '如果达到最大迭代次数，则退出迭代
    If zongi > 300 Then Exit Do
    
Loop Until pp < xe '如果收敛标准的值小于允许误差，则退出迭代
End Sub
Sub duquhuanjing2()
'读取环境中的参数
'读取表格中给如入的3组初始的参数值，存入ps1，ps2，ps3变量中去
Dim p, i, ri
'p参数的个数
p = ActiveSheet.Cells(1000, 9).End(xlUp).Row - 1
ri = p - 1
'遍历表格，索取所有的参数，这里虽然是遍历，但固定死了为3个
For i = 1 To 5
    Select Case i
    Case Is = 1
        Set ps1 = ActiveSheet.Range(ActiveSheet.Cells(2, 9 + i), ActiveSheet.Cells(2 + ri, 9 + i))
    Case Is = 2
        Set ps2 = ActiveSheet.Range(ActiveSheet.Cells(2, 9 + i), ActiveSheet.Cells(2 + ri, 9 + i))
    Case Is = 3
        Set ps3 = ActiveSheet.Range(ActiveSheet.Cells(2, 9 + i), ActiveSheet.Cells(2 + ri, 9 + i))
    '下面两个，nx和n4x是中间结果，被安放在这里指定他俩的存储位置
    Case Is = 4
        Set nx = ActiveSheet.Range(ActiveSheet.Cells(2, 9 + i), ActiveSheet.Cells(2 + ri, 9 + i))
    Case Is = 5
        Set n3x = ActiveSheet.Range(ActiveSheet.Cells(2, 9 + i), ActiveSheet.Cells(2 + ri, 9 + i))
        Set n4x = ActiveSheet.Range(ActiveSheet.Cells(2, 9 + i), ActiveSheet.Cells(2 + ri, 9 + i))
    End Select
Next
'p0 = 目标函数 ，p1=参数1 ，p2=参数2

End Sub

Sub d1_px(ps1, ps2, ps3)
'为3组初始参数进行排序，以便接下来的计算，直接使用其中的最大值与最小值
'fh最大的哪一组,
'fg中间的那一组,
'fl最小的那一组

'第一步，排序
If ps1(1) > ps2(1) Then
    If ps1(1) > ps3(1) Then
    'ps1最大
        Set fh = ps1
        If ps2(1) > ps3(1) Then '判断ps2和ps3
            Set fg = ps2
            Set fl = ps3
        Else
            Set fl = ps2
            Set fg = ps3
        End If
        
    Else
        Set fg = ps1 'ps1是中间的
        If ps2(1) > ps3(1) Then
            Set fh = ps2
            Set fl = ps3
        Else
            Set fl = ps2
            Set fh = ps3
        End If
                
    End If
Else
    If ps1(1) > ps3(1) Then
    'ps1是中间的
        Set fg = ps1
        Set fh = ps2
        Set fl = ps3
    Else
    'ps1是最小的
        Set fl = ps1
        If ps2(1) > ps3(1) Then
            Set fh = ps2
            Set fg = ps3
        Else
            Set fg = ps2
            Set fh = ps3
        End If
    End If
     
End If

'这是无用的步骤用于debug
Dim a, b, c
a = zongi
a = fl(0)
b = fg(0)
c = fh(0)

End Sub

Sub d2_mean(fh, fg, fl)
'计算出fh之外，剩余机组参数中，每个参数的中心点
Dim mx(3), i
For i = 2 To 3
    mx(i) = (fg(i) + fl(i)) / 2
Next

'利用fh和和剩余点的质心，求算反射点
'Dim nx As Range
'Set nx = ActiveSheet.Range(ActiveSheet.Cells(2, 13), ActiveSheet.Cells(4, 13))
For i = 2 To 3
    nx(i) = mx(i) + xa * (mx(i) - fh(i))
Next

'对反射点的大小进行判断，依据不同的情况，求算出新的点或使用质心，作为fh的替代参数
If nx(1) < fl(1) Then
'第5.1步
    'Dim n3x As Range
    'Set n3x = ActiveSheet.Range(ActiveSheet.Cells(2, 14), ActiveSheet.Cells(4, 14))
    For i = 2 To 3
        n3x(i) = mx(i) + xr * (nx(i) - fg(i)) '扩展
    Next
    
    '第6步
    If n3x(1) < nx(1) Then
        For i = 2 To 3
            fh(i) = n3x(i) '直通第8步
        Next
    Else
        For i = 2 To 3
            fh(i) = nx(i) '直通第8步
        Next
    End If
    
Else
    If nx(1) > fg(1) Then
    '第5.3步
        If nx(1) < fh(1) Then
            For i = 2 To 3
                fh(i) = nx(i)
            Next
        End If
        'Dim n4x As Range
        'Set n4x = ActiveSheet.Range(ActiveSheet.Cells(2, 14), ActiveSheet.Cells(4, 14))
        For i = 2 To 3
            n4x(i) = mx(i) + xb * (fh(i) - mx(i))
        Next
        '第7步
        If n4x(1) < fh(1) Then
            'Set fh = n4x '直通第8步，第8步是检验是否满足收敛标准
            For i = 2 To 3
                fh(i) = n4x(i) '直通第8步
            Next
        Else
            For i = 2 To 3
                fh(i) = fh(i) + 0.5 * (fl(i) - fh(i)) '直通第8步
                fg(i) = fg(i) + 0.5 * (fl(i) - fg(i)) '直通第8步
            Next
        End If
    Else
    '第5.2步
        'Set fh = nx '直通第8步
        For i = 2 To 3
            fh(i) = nx(i) '直通第8步
        Next
    End If
End If

End Sub

Function panju(fh, fg, fl)
'第8步，收敛标准值得计算
Dim i, mx(3)
For i = 2 To 3
    mx(i) = (fg(i) + fl(i)) / 2
Next
panju = (( _
(fh(1) - mx(1)) ^ 2 + _
(fg(1) - mx(1)) ^ 2 + _
(fl(1) - mx(1)) ^ 2 _
) / 4) ^ 0.5
End Function

对于方差的计算步骤，我把他放在了Excel表格里来实现，没有通过代码去实现。

Excel表格结构如下图： VBA

自己本地测试通过，可以分非常好的进行回归计。