Use PSO to find minimum in OpenCASCADE
eryar@163.com
Abstract. Starting from OCCT6.8.0 will include one more algorithm for solving global optimization problems. Its development has been triggered by insufficient performance and robustness of existing algorithm of minimization of curve-surface distance in Extrema package. The PSO, Algorithms in this family are stochastic, and this feature can be perceived as opposite to robustness. However, we found it was not only much faster than original deterministic one, but also more robust in complex real-world situations. In particular, it has been able to find solution in situations like tangential or degenerated geometries where deterministic algorithms work poor and require extensive oversampling for robust results. The paper mainly focus on the usage and applications of the PSO algorithm.
Key Words. PSO, Particle Swarm Optimization, Minimization
1.Introduction
粒子群優(yōu)化(Particle Swarm Optimization, PSO)算法是Kennedy和Eberhart受人工生命研究結(jié)果的啟發(fā)、通過(guò)模擬鳥(niǎo)群覓食過(guò)程中的遷徙和群聚行為而提出的一種基于群體智能的全局隨機(jī)搜索算法，自然界中各種生物體均具有一定的群體行為，而人工生命的主要研究領(lǐng)域之一是探索自然界生物的群體行為，從而在計(jì)算機(jī)上構(gòu)建其群體模型。自然界中的鳥(niǎo)群和魚(yú)群的群體行為一直是科學(xué)家的研究興趣，生物學(xué)家Craig Reynolds在1987年提出了一個(gè)非常有影響的鳥(niǎo)群聚集模型，在他的仿真中，每一個(gè)個(gè)體遵循：
(1) 避免與鄰域個(gè)體相沖撞；
(2) 匹配鄰域個(gè)體的速度；
(3) 飛向鳥(niǎo)群中心，且整個(gè)群體飛向目標(biāo)。
仿真中僅利用上面三條簡(jiǎn)單的規(guī)則，就可以非常接近的模擬出鳥(niǎo)群飛行的現(xiàn)象。1995年，美國(guó)社會(huì)心理學(xué)家James Kennedy和電氣工程師Russell Eberhart共同提出了粒子群算法，其基本思想是受對(duì)鳥(niǎo)類(lèi)群體行為進(jìn)行建模與仿真的研究結(jié)果的啟發(fā)。他們的模型和仿真算法主要對(duì)Frank Heppner的模型進(jìn)行了修正，以使粒子飛向解空間并在最好解處降落。Kennedy在他的書(shū)中描述了粒子群算法思想的起源。
粒子群優(yōu)化由于其算法簡(jiǎn)單，易于實(shí)現(xiàn)，無(wú)需梯度信息，參數(shù)少等特點(diǎn)在連續(xù)優(yōu)化問(wèn)題和離散問(wèn)題中都表現(xiàn)出良好的效果，特別是因?yàn)槠涮烊坏膶?shí)數(shù)編碼特點(diǎn)適合處理實(shí)優(yōu)化問(wèn)題。近年來(lái)成為國(guó)際上智能優(yōu)化領(lǐng)域研究的熱點(diǎn)。PSO算法最早應(yīng)用于非線性連續(xù)函數(shù)的優(yōu)化和神經(jīng)元網(wǎng)絡(luò)的訓(xùn)練，后來(lái)也被用于解決約束優(yōu)化問(wèn)題、多目標(biāo)優(yōu)化問(wèn)題，動(dòng)態(tài)優(yōu)化問(wèn)題等。
在OpenCASCADE中很多問(wèn)題可以歸結(jié)為非線性連續(xù)函數(shù)的優(yōu)化問(wèn)題，如求極值的包中計(jì)算曲線和曲面之間的極值點(diǎn)，或曲線與曲線間的極值點(diǎn)等問(wèn)題。本文主要關(guān)注OpenCASCADE中PSO的用法，在理解其用法的基礎(chǔ)上再來(lái)理解其他相關(guān)的應(yīng)用。
2. PSO Usage
為了提高程序的性能及穩(wěn)定性，OpenCASCADE引入了智能化的算法PSO（Particle Swarm Optimization），相關(guān)的類(lèi)為math_PSO，其類(lèi)聲明的代碼如下：
//! In this class implemented variation of Particle Swarm Optimization (PSO) method.
//! A. Ismael F. Vaz, L. N. Vicente
//! "A particle swarm pattern search method for bound constrained global optimization"
//!
//! Algorithm description:
//! Init Section:
//! At start of computation a number of "particles" are placed in the search space.
//! Each particle is assigned a random velocity.
//!
//! Computational loop:
//! The particles are moved in cycle, simulating some "social" behavior, so that new position of
//! a particle on each step depends not only on its velocity and previous path, but also on the
//! position of the best particle in the pool and best obtained position for current particle.
//! The velocity of the particles is decreased on each step, so that convergence is guaranteed.
//!
//! Algorithm output:
//! Best point in param space (position of the best particle) and value of objective function.
//!
//! Pros:
//! One of the fastest algorithms.
//! Work over functions with a lot local extremums.
//! Does not require calculation of derivatives of the functional.
//!
//! Cons:
//! Convergence to global minimum not proved, which is a typical drawback for all stochastic algorithms.
//! The result depends on random number generator.
//!
//! Warning: PSO is effective to walk into optimum surrounding, not to get strict optimum.
//! Run local optimization from pso output point.
//! Warning: In PSO used fixed seed in RNG, so results are reproducible.
class math_PSO
{
public:
/**
* Constructor.
*
* @param theFunc defines the objective function. It should exist during all lifetime of class instance.
* @param theLowBorder defines lower border of search space.
* @param theUppBorder defines upper border of search space.
* @param theSteps defines steps of regular grid, used for particle generation.
This parameter used to define stop condition (TerminalVelocity).
* @param theNbParticles defines number of particles.
* @param theNbIter defines maximum number of iterations.
*/
Standard_EXPORT math_PSO(math_MultipleVarFunction* theFunc,
const math_Vector& theLowBorder,
const math_Vector& theUppBorder,
const math_Vector& theSteps,
const Standard_Integer theNbParticles = 32,
const Standard_Integer theNbIter = 100);
//! Perform computations, particles array is constructed inside of this function.
Standard_EXPORT void Perform(const math_Vector& theSteps,
Standard_Real& theValue,
math_Vector& theOutPnt,
const Standard_Integer theNbIter = 100);
//! Perform computations with given particles array.
Standard_EXPORT void Perform(math_PSOParticlesPool& theParticles,
Standard_Integer theNbParticles,
Standard_Real& theValue,
math_Vector& theOutPnt,
const Standard_Integer theNbIter = 100);
private:
void performPSOWithGivenParticles(math_PSOParticlesPool& theParticles,
Standard_Integer theNbParticles,
Standard_Real& theValue,
math_Vector& theOutPnt,
const Standard_Integer theNbIter = 100);
math_MultipleVarFunction *myFunc;
math_Vector myLowBorder; // Lower border.
math_Vector myUppBorder; // Upper border.
math_Vector mySteps; // steps used in PSO algorithm.
Standard_Integer myN; // Dimension count.
Standard_Integer myNbParticles; // Particles number.
Standard_Integer myNbIter;
};
math_PSO的輸入主要為：求極小值的多元函數(shù)math_MultipleVarFunction，各個(gè)自變量的取值范圍。下面通過(guò)一個(gè)具體的例子來(lái)說(shuō)明如何將數(shù)學(xué)問(wèn)題轉(zhuǎn)換成代碼，利用程序來(lái)求解。在《最優(yōu)化方法》中找到如下例題：
實(shí)現(xiàn)代碼如下所示：
/*
Copyright(C) 2017 Shing Liu(eryar@163.com)
Permission is hereby granted, free of charge, to any person obtaining a copy
of this software and associated documentation files(the "Software"), to deal
in the Software without restriction, including without limitation the rights
to use, copy, modify, merge, publish, distribute, sublicense, and / or sell
copies of the Software, and to permit persons to whom the Software is
furnished to do so, subject to the following conditions :
The above copyright notice and this permission notice shall be included in all
copies or substantial portions of the Software.
THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.IN NO EVENT SHALL THE
AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
SOFTWARE.
*/
// NOTE
// ----
// Tool: Visual Studio 2013 & OpenCASCADE7.1.0
// Date: 2017-04-18 20:52
#include <math_PSO.hxx>
#pragma comment(lib, "TKernel.lib")
#pragma comment(lib, "TKMath.lib")
// define function:
// f(x) = x1 - x2 + 2x1^2 + 2x1x2 + x2^2
class TestFunction : public math_MultipleVarFunction
{
public:
virtual Standard_Integer NbVariables() const
{
return 2;
}
virtual Standard_Boolean Value(const math_Vector& X, Standard_Real& F)
{
F = X(1) - X(2) + 2.0 * X(1) * X(1) + 2.0 * X(1) * X(2) + X(2) * X(2);
return Standard_True;
}
};
void test()
{
TestFunction aFunction;
math_Vector aLowerBorder(1, aFunction.NbVariables());
math_Vector aUpperBorder(1, aFunction.NbVariables());
math_Vector aSteps(1, aFunction.NbVariables());
aLowerBorder(1) = -10.0;
aLowerBorder(2) = -10.0;
aUpperBorder(1) = 10.0;
aUpperBorder(2) = 10.0;
aSteps(1) = 0.1;
aSteps(2) = 0.1;
Standard_Real aValue = 0.0;
math_Vector aOutput(1, aFunction.NbVariables());
math_PSO aPso(&aFunction, aLowerBorder, aUpperBorder, aSteps);
aPso.Perform(aSteps, aValue, aOutput);
std::cout << "Minimum value: " << aValue << " at (" << aOutput(1) << ", " << aOutput(2) << ")" << std::endl;
}
int main(int argc, char* argv[])
{
test();
return 0;
}
計(jì)算結(jié)果為最小值-1.25, 在x1=-1，x2=1.50013處取得，如下圖所示：
3.Applications
計(jì)算曲線和曲面之間的極值，也是一個(gè)計(jì)算多元函數(shù)的極值問(wèn)題。對(duì)應(yīng)的類(lèi)為Extrema_GenExtCS。
4.Conclusion
常見(jiàn)的優(yōu)化算法如共軛梯度法、Netwon-Raphson法（math_NewtonMinimum），F(xiàn)letcher-Reeves法（math_FRPR）、Broyden-Fletcher-Glodfarb-Shanno法（math_BFGS）等都是局部?jī)?yōu)化方法，所有這些局部?jī)?yōu)化方法都是針對(duì)無(wú)約束優(yōu)化問(wèn)題提出的，而且對(duì)目標(biāo)函數(shù)均有一定的解析性要求，如Newton-Raphson要求目標(biāo)函數(shù)連續(xù)可微，同時(shí)要求其一階導(dǎo)數(shù)連續(xù)。OpenCASCADE中的math_NewtonMinimu中還要求多元函數(shù)具有Hessian矩陣，即二階導(dǎo)數(shù)連續(xù)。
隨著現(xiàn)科學(xué)技術(shù)的不斷發(fā)展和多學(xué)科的相互交叉與滲透，很多實(shí)際工程優(yōu)化問(wèn)題不僅需要大量的復(fù)雜科學(xué)計(jì)算，而且對(duì)算法的實(shí)時(shí)性要求也特別高，這些復(fù)雜優(yōu)化問(wèn)題通常具有以下特點(diǎn)：
1) 優(yōu)化問(wèn)題的對(duì)象涉及很多因素，導(dǎo)致優(yōu)化問(wèn)題的目標(biāo)函數(shù)的自變量維數(shù)很多，通常達(dá)到數(shù)百維甚至上萬(wàn)維，使得求解問(wèn)題的計(jì)算量大大增加；
2) 優(yōu)化問(wèn)題本身的復(fù)雜性導(dǎo)致目標(biāo)函數(shù)是非線性，同時(shí)由于有些目標(biāo)函數(shù)具有不可導(dǎo)，不連續(xù)、極端情況下甚至函數(shù)本身不能解析的表達(dá)；如曲線或曲面退化導(dǎo)致的尖點(diǎn)或退化點(diǎn)的不連續(xù)情況；
3) 目標(biāo)函數(shù)在定義域內(nèi)具有多個(gè)甚至無(wú)數(shù)個(gè)極值點(diǎn)，函數(shù)的解空間形狀復(fù)雜。
當(dāng)以上三個(gè)因素中一個(gè)或幾個(gè)出現(xiàn)在優(yōu)化問(wèn)題中時(shí)，將會(huì)極大地增加優(yōu)化問(wèn)題求解的困難程度，單純地基于解析確定性的優(yōu)化方法很難奏效，因此必須結(jié)合其他方法來(lái)解決這些問(wèn)題。
PSO算法簡(jiǎn)單，易于實(shí)現(xiàn)且不需要目標(biāo)函數(shù)的梯度（一階可導(dǎo)）和Hassian矩陣（二階可導(dǎo)），速度快，性能高。但是PSO也有不足之處，這是許多智能算法的共性，即只能找到滿足條件的解，這個(gè)解不能確定為精確解。
5.References
1. aml. Application of stohastic algorithms in extrema. https://dev.opencascade.org/index.php?q=node/988
2. 何堅(jiān)勇. 最優(yōu)化方法. 清華大學(xué)出版社. 2007
3. 沈顯君. 自適應(yīng)粒子群優(yōu)化算法及其應(yīng)用. 清華大學(xué)出版社. 2015
4. 汪定偉, 王俊偉, 王洪峰, 張瑞友, 郭哲. 智能優(yōu)化方法. 高等教育出版社. 2007