精华内容
下载资源
问答
  • 2021-04-07 11:29:02

    ORB-SLAM2学习笔记——局部BA优化

    1、理论部分(待更新)

    2、代码详解

    void Optimizer::LocalBundleAdjustment(KeyFrame *pKF, bool* pbStopFlag, Map* pMap)
    {
        // 该优化函数用于LocalMapping线程的局部BA优化    
        //step 1 : 得到所有局部关键帧
        //step 2 : 得到所有地图点
        //step 3 : 把关键帧和能看到的地图点对应到一起
        //step 4 : 进行优化 
    
        // Local KeyFrames: First Breath Search from Current Keyframe
        //局部 关键帧: 当前帧的第一次轨迹搜索
        //定义局部关键帧的列表
        list<KeyFrame*> lLocalKeyFrames;
        //关键帧放入列表
        lLocalKeyFrames.push_back(pKF);
        //索引导入到匹配好的局部优化关键帧
        pKF->mnBALocalForKF = pKF->mnId;
        //把匹配上的关键帧的顺序排好进行保存
        const vector<KeyFrame*> vNeighKFs = pKF->GetVectorCovisibleKeyFrames();
        //遍历关键帧保存到pKFi中,然后把相应的索引也赋值过来
        for(int i=0, iend=vNeighKFs.size(); i<iend; i++)
        {
            KeyFrame* pKFi = vNeighKFs[i];
            pKFi->mnBALocalForKF = pKF->mnId;
            if(!pKFi->isBad())
                lLocalKeyFrames.push_back(pKFi);
        }
    
        // Local MapPoints seen in Local KeyFrames
        //被局部关键帧看见的局部地图点
        list<MapPoint*> lLocalMapPoints;
        //遍历所有地图点
        for(list<KeyFrame*>::iterator lit=lLocalKeyFrames.begin() , lend=lLocalKeyFrames.end(); lit!=lend; lit++)
        {
            //取出匹配上的地图点
            vector<MapPoint*> vpMPs = (*lit)->GetMapPointMatches();
            //遍历所有匹配上的地图点剔除坏点,保存好点及其索引
            for(vector<MapPoint*>::iterator vit=vpMPs.begin(), vend=vpMPs.end(); vit!=vend; vit++)
            {
                MapPoint* pMP = *vit;
                if(pMP)
                    if(!pMP->isBad())
                        if(pMP->mnBALocalForKF!=pKF->mnId)
                        {
                            lLocalMapPoints.push_back(pMP);
                            pMP->mnBALocalForKF=pKF->mnId;// 防止重复添加
                        }
            }
        }
    
        // Fixed Keyframes. Keyframes that see Local MapPoints but that are not Local Keyframes
        //得到能被局部MapPoints观测到,但不属于局部关键帧的关键帧,这些关键帧在局部BA优化时不优化
        //筛选出能看见对应局部地图点的局部关键帧
        list<KeyFrame*> lFixedCameras;
        //遍历所有局部地图点
        for(list<MapPoint*>::iterator lit=lLocalMapPoints.begin(), lend=lLocalMapPoints.end(); lit!=lend; lit++)
        {   //得到对应的观测
            map<KeyFrame*,size_t> observations = (*lit)->GetObservations();
            //遍历所有的观测,得到对应的关键帧,剔除坏帧,保存对应索引
            for(map<KeyFrame*,size_t>::iterator mit=observations.begin(), mend=observations.end(); mit!=mend; mit++)
            {
                KeyFrame* pKFi = mit->first;
    
                if(pKFi->mnBALocalForKF!=pKF->mnId && pKFi->mnBAFixedForKF!=pKF->mnId)
                {                
                    pKFi->mnBAFixedForKF=pKF->mnId;
                    if(!pKFi->isBad())
                        lFixedCameras.push_back(pKFi);
                }
            }
        }
    
        // Setup optimizer
        //创建优化器
        g2o::SparseOptimizer optimizer;
        g2o::BlockSolver_6_3::LinearSolverType * linearSolver;
    
        linearSolver = new g2o::LinearSolverEigen<g2o::BlockSolver_6_3::PoseMatrixType>();
    
        g2o::BlockSolver_6_3 * solver_ptr = new g2o::BlockSolver_6_3(linearSolver);
        //选用LM算法进行求解线性方程
        g2o::OptimizationAlgorithmLevenberg* solver = new g2o::OptimizationAlgorithmLevenberg(solver_ptr);
        optimizer.setAlgorithm(solver);
        //停止优化标志
        if(pbStopFlag)
            optimizer.setForceStopFlag(pbStopFlag);
    
        unsigned long maxKFid = 0;
    
        // Set Local KeyFrame vertices
        //创建局部关键帧的顶点
        for(list<KeyFrame*>::iterator lit=lLocalKeyFrames.begin(), lend=lLocalKeyFrames.end(); lit!=lend; lit++)
        {
            KeyFrame* pKFi = *lit;
    
            g2o::VertexSE3Expmap * vSE3 = new g2o::VertexSE3Expmap();
            //设置待估计值 
            vSE3->setEstimate(Converter::toSE3Quat(pKFi->GetPose()));
            //设置索引
            vSE3->setId(pKFi->mnId);
            //固定第一帧防止整体乱动
            vSE3->setFixed(pKFi->mnId==0);
            //添加顶点到优化器
            optimizer.addVertex(vSE3);
            
            if(pKFi->mnId>maxKFid)
                maxKFid=pKFi->mnId;
        }
    
        // Set Fixed KeyFrame vertices
        //建立要固定的关键帧顶点 
        for(list<KeyFrame*>::iterator lit=lFixedCameras.begin(), lend=lFixedCameras.end(); lit!=lend; lit++)
        {
            KeyFrame* pKFi = *lit;
            g2o::VertexSE3Expmap * vSE3 = new g2o::VertexSE3Expmap();
            vSE3->setEstimate(Converter::toSE3Quat(pKFi->GetPose()));
            vSE3->setId(pKFi->mnId);
            //不进行优化
            vSE3->setFixed(true);
            optimizer.addVertex(vSE3);
            if(pKFi->mnId>maxKFid)
                maxKFid=pKFi->mnId;
        }
    
        // Set MapPoint vertices
        ///添加3D顶点
        const int nExpectedSize = (lLocalKeyFrames.size()+lFixedCameras.size()) * lLocalMapPoints.size();
    
        vector<g2o::EdgeSE3ProjectXYZ*> vpEdgesMono;
        vpEdgesMono.reserve(nExpectedSize);
    
        vector<KeyFrame*> vpEdgeKFMono;
        vpEdgeKFMono.reserve(nExpectedSize);
    
        vector<MapPoint*> vpMapPointEdgeMono;
        vpMapPointEdgeMono.reserve(nExpectedSize);
    
        vector<g2o::EdgeStereoSE3ProjectXYZ*> vpEdgesStereo;
        vpEdgesStereo.reserve(nExpectedSize);
    
        vector<KeyFrame*> vpEdgeKFStereo;
        vpEdgeKFStereo.reserve(nExpectedSize);
    
        vector<MapPoint*> vpMapPointEdgeStereo;
        vpMapPointEdgeStereo.reserve(nExpectedSize);
        //鲁棒核系数
        const float thHuberMono = sqrt(5.991);
        const float thHuberStereo = sqrt(7.815);
        //添加局部地图点到顶点
        for(list<MapPoint*>::iterator lit=lLocalMapPoints.begin(), lend=lLocalMapPoints.end(); lit!=lend; lit++)
        {
            MapPoint* pMP = *lit;
            g2o::VertexSBAPointXYZ* vPoint = new g2o::VertexSBAPointXYZ();
            vPoint->setEstimate(Converter::toVector3d(pMP->GetWorldPos()));
            int id = pMP->mnId+maxKFid+1;
            vPoint->setId(id);
            //进行边缘化
            vPoint->setMarginalized(true);
            optimizer.addVertex(vPoint);
    
            const map<KeyFrame*,size_t> observations = pMP->GetObservations();
    
            //Set edges
            //创建边
    
            //遍历所有观测
            for(map<KeyFrame*,size_t>::const_iterator mit=observations.begin(), mend=observations.end(); mit!=mend; mit++)
            {   
                //取出观测对应关键帧
                KeyFrame* pKFi = mit->first;
                //剔除坏点
                if(!pKFi->isBad())
                {                
                    const cv::KeyPoint &kpUn = pKFi->mvKeysUn[mit->second];
    
                    // Monocular observation
                    //单目观测
                    if(pKFi->mvuRight[mit->second]<0)
                    {
                        //观测值
                        Eigen::Matrix<double,2,1> obs;
                        obs << kpUn.pt.x, kpUn.pt.y;
    
                        g2o::EdgeSE3ProjectXYZ* e = new g2o::EdgeSE3ProjectXYZ();
                        //优化的顶点 0 地图点
                        e->setVertex(0, dynamic_cast<g2o::OptimizableGraph::Vertex*>(optimizer.vertex(id)));
                        //优化的顶点 1 位姿
                        e->setVertex(1, dynamic_cast<g2o::OptimizableGraph::Vertex*>(optimizer.vertex(pKFi->mnId)));
                        //观测值
                        e->setMeasurement(obs);
                        //取出对应点的权重
                        const float &invSigma2 = pKFi->mvInvLevelSigma2[kpUn.octave];
                        //建立信息矩阵
                        e->setInformation(Eigen::Matrix2d::Identity()*invSigma2);
    
                        //建立鲁棒核函数
                        g2o::RobustKernelHuber* rk = new g2o::RobustKernelHuber;
                        e->setRobustKernel(rk);
                        rk->setDelta(thHuberMono);
                        //相机内参
                        e->fx = pKFi->fx;
                        e->fy = pKFi->fy;
                        e->cx = pKFi->cx;
                        e->cy = pKFi->cy;
                        //添加边到优化
                        optimizer.addEdge(e);
                        vpEdgesMono.push_back(e);
                        vpEdgeKFMono.push_back(pKFi);
                        vpMapPointEdgeMono.push_back(pMP);
                    }
                    else // Stereo observation 双目
                    {
                        Eigen::Matrix<double,3,1> obs;
                        const float kp_ur = pKFi->mvuRight[mit->second];
                        obs << kpUn.pt.x, kpUn.pt.y, kp_ur;
    
                        g2o::EdgeStereoSE3ProjectXYZ* e = new g2o::EdgeStereoSE3ProjectXYZ();
    
                        e->setVertex(0, dynamic_cast<g2o::OptimizableGraph::Vertex*>(optimizer.vertex(id)));
                        e->setVertex(1, dynamic_cast<g2o::OptimizableGraph::Vertex*>(optimizer.vertex(pKFi->mnId)));
                        e->setMeasurement(obs);
                        const float &invSigma2 = pKFi->mvInvLevelSigma2[kpUn.octave];
                        Eigen::Matrix3d Info = Eigen::Matrix3d::Identity()*invSigma2;
                        e->setInformation(Info);
    
                        g2o::RobustKernelHuber* rk = new g2o::RobustKernelHuber;
                        e->setRobustKernel(rk);
                        rk->setDelta(thHuberStereo);
    
                        e->fx = pKFi->fx;
                        e->fy = pKFi->fy;
                        e->cx = pKFi->cx;
                        e->cy = pKFi->cy;
                        e->bf = pKFi->mbf;
    
                        optimizer.addEdge(e);
                        vpEdgesStereo.push_back(e);
                        vpEdgeKFStereo.push_back(pKFi);
                        vpMapPointEdgeStereo.push_back(pMP);
                    }
                }
            }
        }
        //
        if(pbStopFlag)
            if(*pbStopFlag)
                return;
        //初始化优化器,开始优化
        optimizer.initializeOptimization();
        optimizer.optimize(5);
    
        bool bDoMore= true;
        //终止优化标志
        if(pbStopFlag)
            if(*pbStopFlag)
                bDoMore = false;
        //
        if(bDoMore)
        {
    
        // Check inlier observations
        //检测outlier,并设置下次不优化
        for(size_t i=0, iend=vpEdgesMono.size(); i<iend;i++)
        {
            g2o::EdgeSE3ProjectXYZ* e = vpEdgesMono[i];
            MapPoint* pMP = vpMapPointEdgeMono[i];
            //如果是坏帧跳过
            if(pMP->isBad())
                continue;
            //误差大于阈值或深度没有测量到
            if(e->chi2()>5.991 || !e->isDepthPositive())
            {
                //设置为标志 1 不进行优化
                e->setLevel(1);
            }
            //鲁棒核函数设置为 0 不使用核函数
            e->setRobustKernel(0);
        }
    
        for(size_t i=0, iend=vpEdgesStereo.size(); i<iend;i++)
        {
            g2o::EdgeStereoSE3ProjectXYZ* e = vpEdgesStereo[i];
            MapPoint* pMP = vpMapPointEdgeStereo[i];
    
            if(pMP->isBad())
                continue;
    
            if(e->chi2()>7.815 || !e->isDepthPositive())
            {
                e->setLevel(1);
            }
    
            e->setRobustKernel(0);
        }
    
        // Optimize again without the outliers
        //不带外点进行再次优化
        optimizer.initializeOptimization(0);
        optimizer.optimize(10);
    
        }
        //
        vector<pair<KeyFrame*,MapPoint*> > vToErase;
        vToErase.reserve(vpEdgesMono.size()+vpEdgesStereo.size());
    
        // Check inlier observations  
        //降低误差阈值,再次进行筛选     
        for(size_t i=0, iend=vpEdgesMono.size(); i<iend;i++)
        {
            g2o::EdgeSE3ProjectXYZ* e = vpEdgesMono[i];
            MapPoint* pMP = vpMapPointEdgeMono[i];
    
            if(pMP->isBad())
                continue;
            // 基于卡方检验计算出的阈值(假设测量有一个像素的偏差)
            if(e->chi2()>5.991 || !e->isDepthPositive())
            {
                KeyFrame* pKFi = vpEdgeKFMono[i];
                //要去除关键帧放在一起
                vToErase.push_back(make_pair(pKFi,pMP));
            }
        }
        //双目的
        for(size_t i=0, iend=vpEdgesStereo.size(); i<iend;i++)
        {
            g2o::EdgeStereoSE3ProjectXYZ* e = vpEdgesStereo[i];
            MapPoint* pMP = vpMapPointEdgeStereo[i];
    
            if(pMP->isBad())
                continue;
    
            if(e->chi2()>7.815 || !e->isDepthPositive())
            {
                KeyFrame* pKFi = vpEdgeKFStereo[i];
                vToErase.push_back(make_pair(pKFi,pMP));
            }
        }
    
        // Get Map Mutex
        //
        unique_lock<mutex> lock(pMap->mMutexMapUpdate);
        //删除误差比较大的关键帧和关键点
        if(!vToErase.empty())
        {
            for(size_t i=0;i<vToErase.size();i++)
            {
                KeyFrame* pKFi = vToErase[i].first;
                MapPoint* pMPi = vToErase[i].second;
                pKFi->EraseMapPointMatch(pMPi);
                pMPi->EraseObservation(pKFi);
            }
        }
    
        // Recover optimized data
        //恢复出优化后的数据
        //Keyframes
        //恢复出关键帧的信息
        for(list<KeyFrame*>::iterator lit=lLocalKeyFrames.begin(), lend=lLocalKeyFrames.end(); lit!=lend; lit++)
        {
            KeyFrame* pKF = *lit;
            g2o::VertexSE3Expmap* vSE3 = static_cast<g2o::VertexSE3Expmap*>(optimizer.vertex(pKF->mnId));
            g2o::SE3Quat SE3quat = vSE3->estimate();
            pKF->SetPose(Converter::toCvMat(SE3quat));
        }
    
        //Points
        //恢复出关键点的信息
        for(list<MapPoint*>::iterator lit=lLocalMapPoints.begin(), lend=lLocalMapPoints.end(); lit!=lend; lit++)
        {
            MapPoint* pMP = *lit;
            g2o::VertexSBAPointXYZ* vPoint = static_cast<g2o::VertexSBAPointXYZ*>(optimizer.vertex(pMP->mnId+maxKFid+1));
            pMP->SetWorldPos(Converter::toCvMat(vPoint->estimate()));
            pMP->UpdateNormalAndDepth();
        }
    }
    
    
    更多相关内容
  • 通过约束局部优化算子的参数,减少局部优化的计算量,使整体算法的复杂度基本蚁群算法大致相当.从最终的实验数据可以得出,所提出的算法在较少迭代次数的情况下可以得出较高的精度,在收敛速度收敛精度之间实现较好...
  • 局部优化模型中, 设计区域对比度描述子和区域特征描述子对多级分割的超像素点进行描述, 预测每一个区域的显著性值。最后, 利用线性拟合的方法将两种模型中产生的显著图进行融合, 得到最终的显著图。对4个数据集...
  • 转运作为主要的站内作业活动,对其开展调度研究,能够有效缩短不同运输方式之间的转运周期,从而保证铁路集装箱多式联运的整体运作效率.为了快速制定堆场作业计划,现有铁路集装箱站场转运作业均要求班列间的集装箱...
  • 在JavaScript中,数据的存储位置对代码的整体性能有着重要的影响。有四种数据访问类型:直接量,局部变量,数组项,对象成员
  • 然后,为了很好地平衡探索开采能力以提升整体优化性能,对算法前、后半搜索阶段分别采用单维操作和全维操作形成ODGWO;最后,将ODGWO用于高维函数和模糊C均值(FCM)聚类优化.实验结果表明,在许多高维Benchmark函数(30...
  • 针对局部立体匹配算法在弱纹理区域匹配准确度低的问题,提出一种在代价初始化和代价聚合两阶段自适应优化的立体匹配算法。首先,在代价初始化阶段,通过双阈值线性约束条件构造像素点的十字支撑窗口,并依据十字支撑最短...
  • 现阶段NB-IoT网络已完成初步建设,整体覆盖较好,但局部覆盖空洞较多,室内疑难场景深度覆盖不足,干扰问题突出,重叠覆盖严重。目前NB-IoT行业应用成熟度低,行业应用属于起步阶段,各厂商终端规范不一,导致网络...
  • 首先,采用半均匀初始化种群,使种群以整体均匀、局部随机的方式分布;其次,引入动态分裂算子,对满足分裂条件的粒子执行分裂操作,增加种群多样性,避免粒子陷入局部最优;最后,采用指数衰减的惯性权重,平衡粒子...
  • , 首先教你SQL整体优化、快速优化实施、如何读懂执行计划、如何左右执行计划这四大必杀招。整这些干嘛呢?答案是,传授一个先整体后局部的宏观解决思路,走进“道”的世界。, 接下来带领大家飞翔在“术”的天空。教...
  • 上一次记录了解决过度绘制的...对于顺畅度的分析,首先要知道一个整体情况,是局部,还是全局,这样在优化上才能有方向。如果是局部问题,那就需要仔细分析出具体的相关操作,如果是大体上的问题,那在思考的时候,就需
  • 利用对立学习机制增加初始群体的多样性,并引入和声搜索以增强局部寻优能力, 从而提升整体寻优性能.以典型混沌系统为例进行了未知参数估计的数值仿真,结果验证了所提出混合生物地理优化方法的有效性和鲁棒性.
  • 针对混沌蚂蚁群优化算法(CASO) 容易陷入局部极值和精度低的缺陷, 从认知学角度进行分析, 将创造性思维(CT) 引入CASO 算法, 提出了一种带创造性思维的混沌蚂蚁群优化算法(CTCASO). 基于CT 过程的“四阶段”模型, 构建...
  • 该算法的基本思路是边构造边调整路径,在调整中采用了独创的逆向调整方法,避免算法陷入局部优化陷阱。理论分析和大量实验结果表明,该算法不仅时间复杂度和空间复杂度低,寻优能力也相当强,其综合性能超过目前的一些...
  • 从递归看动态规划前情提要最优问题的本质目标函数最优组合的求解策略枚举递归斐波那契数列问题描述示例解题思路迭代递归...调试问题暴力递归的优化——剪枝从整个搜索策略上来调整从解空间图的角度做调整总结升华...

    前情提要

    在上一篇文章(从贪心算法开始认识动态规划——硬币找零问题)里,我们已经学习了贪心算法的思想,并且发现贪心算法是一种使用局部最优思想解题的算法,即从问题的某一个初始解出发逐步逼近给定的目标,以尽可能快的速度去求得更好的解,当达到算法中的某一步不能再继续前进时,算法停止。

    通过硬币找零问题,我们也发现了贪心算法本身的局限性:不能保证求得的最后解是最佳的

    因此,我们需要从整体最优的角度来解决算法问题。

    最优问题的本质

    所谓最优化问题,就是指在某些约束条件下,决定可选择的变量应该取何值,使所选定的目标函数达到最优的问题。

    从数学意义上说,最优化方法是一种求极值的方法,即在一组约束为等式或不等式的条件下,使系统的目标函数达到极值,即最大值或最小值。在数学里一切都是函数,现在我们先把这个问题用函数形式来表示。

    目标函数

    假定给出 y 元硬币,硬币面额是 5 元和 3 元,求出需要的最少硬币数量。所谓的最少硬币数量就是 5 元硬币和 3 元硬币的总数,假定 5 元硬币数量为 x 0 x_0 x0​,3 元硬币数量为 x 1 x_1 x1​,那么用函数表示就是:

    • f ( x 0 ​ , x 1 ​ ) = x 0 ​ + x 1 f(x_0​,x_1​)=x_0​+x_1 f(x0,x1)=x0+x1

    这就是所谓的“目标函数”

    但是这个函数现在是没有任何限制的,我们希望对此进行约束,使得 5 元硬币和 3 元硬币的面值综合为 y。为此我们需要给出一个约束:

    • 5 x 0 ​ + 3 x 1 ​ = y 5x_0​+3x_1​=y 5x0+3x1=y

    这个时候我们的问题就变成了,当满足这个约束条件的时候,求解函数中的变量 x 0 x_0 x0​ 和 x 1 x_1 x1​,使得目标函数 f ( x 0 ​ , x 1 ​ ) f(x_0​,x_1​) f(x0,x1) 的取值最小。如果用数学的描述方法来说的话,就是下面这样:

    • a r g m i n ( x 0 ​ , x 1 ​ ) ∈ S ​ ( x 0 ​ + x 1 ​ ) argmin(x_0​,x_1​)∈S​(x_0​+x_1​) argmin(x0,x1)S(x0+x1)

    这个就是我们常见的 argmin 表示方式。它的意思是:当 ( x 0 ​ , x 1 ​ ) (x_0​,x_1​) (x0,x1) 属于 S 这个集合的时候,希望知道 x 0 ​ + x 1 x_0​+x_1 x0+x1​ 的最小值是多少。其中 S 集合的条件就是上面的约束。

    所以最优化问题在我们生活中是非常普遍的,只不过大多数问题可能都像硬币找零问题这样看起来普普通通,概念其实是不难理解的。

    回到硬币找零这个问题上。由于 ( x 0 ​ , x 1 ​ ) (x_0​,x_1​) (x0,x1) 都是离散的值,因此所有满足上述约束的 ( x 0 ​ , x 1 ​ ) (x_0​,x_1​) (x0,x1) 组合,就是我们最终所求的集合。

    而这个最优化问题的本质就是:

    • 从所有满足条件的组合 ( x 0 ​ , x 1 ​ ) (x_0​,x_1​) (x0,x1) 中找出一个组合,使得 x 0 ​ + x 1 x_0​+x_1 x0+x1​ 的值最小。

    所以,你会发现在这种离散型的最优化问题中,本质就是从所有满足条件的组合(能够凑出 y 元的组合)中选择出使得我们的目标函数(所有硬币数量之和)最小的那个组合。而这个所谓满足条件的组合不就是 argmin 公式中的那个集合 S 吗?

    因此,这种离散型的最优化问题就是去所有满足条件的组合里找出最优解的组合。

    换句话说,局部最优就是在一定条件下的最优解,而整体最优就是我们真正希望得到的最优解。

    那么,怎么去寻找最优解呢?

    最优组合的求解策略

    枚举

    如果想得到最优组合,那么最简单直接的方法肯定就是枚举。

    枚举就是直接求出所有满足条件的组合,然后看看这些组合是否能得到最大值或者最小值。

    在硬币找零问题中,假设现在需要给出 25 元的硬币,有两种组合,分别是 (5, 0) 和 (2, 5),也就是 5 个 5 元硬币,或者 2 个 5 元硬币加上 5 个 3 元硬币,那么硬币数量最小的组合肯定就是 (5, 0)。

    所以最简单的方法就是找出所有满足条件的组合,也就是上面两个组合,然后去看这些组合中的最优解。枚举本身很简单,就是把所有组合都遍历一遍即可。

    可现在问题就是,如何得到这些组合呢?

    这就需要我们通过一些策略来生成所有满足条件的组合。而递归正是得到这些组合的方法。

    递归

    其实最优化问题使用递归来处理是非常清晰的,递归是搜索组合的一种非常直观的思路。

    在动态规划中,几乎所有问题都需要被描述成递归的形式。

    斐波那契数列

    严格来说,斐波那契数列问题不是最优化问题,但它能很好地展示递归的概念。

    先来看一下斐波那契数列的问题描述。

    问题描述

    斐波那契数通常用 F(n) 表示,形成的序列称为斐波那契数列。该数列由 0 和 1 开始,后面的每一项数字都是前面两项数字的和:
    在这里插入图片描述

    示例

    示例 1:
    
    输入:2
    输出:1
    解释:F(2) = F(1) + F(0) = 1 + 0 = 1。
    
    示例 2:
    
    输入:3
    输出:2
    解释:F(3) = F(2) + F(1) = 1 + 1 = 2.
    

    解题思路

    很多人在解算法面试问题的时候有一种倾向性,那就是使用迭代而非递归来求解问题。

    先不说这样的倾向性正确与否,那么我们就按照这个偏好来解一下(即斐波那契数列的循环解法)。

    迭代
    def fibonacci(n):
        resolution = [0, 1] # 存放结果
        if (n < 2):
            return resolution[n]
    
        # 分别代表第item个值,第item个值的前两个值,第item的值
        item, num1, num2, fib = 1, 0, 1, 0
        while (item < n):
            fib = num1 + num2
            num1 = num2
            num2 = fib
            item += 1
        return fib
    
    
    def main():
        result = fibonacci(4)
        print(result)
    
    if __name__ == "__main__":
        main()
    

    这样的解法固然没错,但是它几乎脱离了题设的数学表达形式。

    在这道题目中,出题者“刻意”地写出了求解斐波那契数列的函数表达式,这其中有没有什么别的含义或原因呢?

    当然有了,这个函数表达式很好地反应出了计算机科学中常见的算法形式:递归。 下面,就让我们来看看斐波那契数列与递归之间的关系。

    递归

    事实上,斐波那契数列的数学形式就是递归的,通过代码你应该能很清楚地看出这一点:

    def Fibonacci(n):
        if (n == 0 or n == 1):
            return n
        if (n > 1):
            return Fibonacci(n - 1) + Fibonacci(n - 2)
    
    def main():
        result = Fibonacci(4)
        print(result)
    
    if __name__ == "__main__":
        main()
    

    递归形式的求解几乎就是简单的把题设中的函数表达式照搬过来,因此从数学意义上讲,递归更直观,且易于理解。

    深入理解递归

    堆栈与递归的状态存储

    在计算机中,实现递归必须建立在堆栈的基础上,这是因为每次递归调用的时候我们都需要把当前函数调用中的局部变量保存在某个特定的地方,等到函数返回的时候再把这些局部变量取出来。

    而用于保存这些局部变量的地方也就是堆栈了。因此,递归可以不断保存当前求解状态并进入下一层次的求解,并在得到后续阶段的解之后,将当前求解状态恢复并与后续求解结果进行合并。

    在硬币找零问题中,我们可以放心的在函数中用循环不断遍历,找出当前面值硬币的可能数量。而无需用其它方法来存储当前或之前的数据。得益于递归,我们通过堆栈实现了状态存储,这样的代码看起来简单、清晰明了。

    递归与回溯

    在求解最优化问题的时候,我们经常会用到回溯这个策略。

    上一篇文章里,我们已经提到过回溯的思想。在硬币找零这个问题里,具体说就是如果遇到已经无法求解的组合,那么我们就往回退一步,修改上一个面值的硬币数量,然后再尝试新的组合。

    递归这种形式,正是赋予了回溯这种可以回退一步的能力:它通过堆栈保存了上一步的当前状态。

    因此,如果想要用回溯的策略来解决问题,那么递归应该是你的首选方法。所以说,回溯在最优化问题中有多么重要,递归也就有多么重要。

    树形结构与深度优先搜索

    为了理解递归,这里用树来描述递归的求解过程:
    在这里插入图片描述

    你可以从中看到形象的递归求解过程:

    • 每个节点的 /(斜线)左边表示当前节点使用的硬币面值,右边表示使用面值后的余额。
    • 蓝色节点就表示我们目前得到的解。

    递归的过程的确就是一个树形结构,而递归也就是一个深度优先搜索的过程,先找到下一步的解,然后再回退,如此往复。

    所以我们可以这样理解递归:作为一个算法解决方案,它采用了深度优先搜索的策略,去搜索所有可能的组合,并得到最优解的最优化问题。

    如果在每个节点上加上当前这个节点求得的组合结果,就可以用递归树表示求解的组合空间

    在这里插入图片描述

    从上图中我们可以发现,每个节点都存储了一个当前求解过程中的组合,和后续节点的组合合并到一起形成完整的答案。

    而真正的求解组合,就是把所有余额为 0 的组合拿出来,经过去重之后得到的结果。

    所以,你可以看到求解的组合就蕴含在这个递归的树形结算的节点空间中,这也就是为什么递归策略是行之有效的:我们可以通过穷举法从所有的解中得到最优解!

    暴力递归的问题

    虽然递归能解决问题,但是这样的穷举法效率实在低下,不仅如此,这样的代码可读性低且调试困难。

    性能问题

    暴力递归的最主要的特点就是穷举(暴力,你说是不是)。

    如果我们只使用朴素的递归思路解题,就需要通过递归来暴力穷举出所有的组合,而且我们穷举的不只是组合,还是所有可能得到目标组合的组成路径

    这个在上面的图中我们可以看到,同样是求解 (2, 5) 这个组合,图中有多少种路径?这还只是 25 元和两种面值的情况。如果求解的金额和面值数量增加,那么我们可以看到这个树会以非常难以置信的方式增长,那么带来的性能问题就是灾难性的。如果你仔细观察一下,就会发现这个树会随着总额的增加呈现指数形式的增长。对于这种事情,我们难以接受。

    因此,递归只是让问题可以求解,但是如果数据规模过大的时候暴力递归会引发极大的性能问题

    可读性与调试问题

    虽然递归在数学意义上非常直观,但是如果问题过于复杂,一般是无法直接画出上面我画的那棵求解树的。

    画求解树的时候,我们可以想出我们的求解过程是怎么进行的,但如果求解树的分支极多,那么很多人就很难继续在脑海中模拟出整个求解过程了。

    因此,一旦程序出现 bug,当你想尝试去调试的时候,就会发现这样的代码几乎没有调试的可能性。这种问题在数据规模很大的情况下尤为明显。

    暴力递归的优化——剪枝

    从前面的图中看到,这棵树中有很多分支是完全相同的,从理论上讲最终只有两个组合,但是这棵树到达同一种组合的路径却非常多,所以优化递归的思路其实就是如何减少搜索的分支数量。

    分支数量减少了,递归效率也就高了,这就是所谓的剪枝优化。

    从整个搜索策略上来调整

    第一种思路是仿照贪心算法,从整个搜索策略上来调整。也就是说,你要考虑这个问题的性质,即面值大的硬币用得足够多,那么这个组合的硬币总数肯定就最小。

    所以在每一次递归时,我们不应该暴力地搜索所有的面值,而应该从面值最大的硬币着手,不断尝试大面值硬币的最大情况。如果无法满足条件再减少一个,再递归搜索。最后的代码就跟上一篇文章中的回溯代码一样,即通过贪心这种思路结合递归实现一种组合搜索。

    殊途同归!从递归的角度重新解释了这个算法问题,而且代码实现也是一样的。

    从解空间图的角度做调整

    在解空间的图中,只要是余额相同的情况下,后面的搜索路径是完全一致的:

    在这里插入图片描述
    在图中圈出的两个部分就是重复的搜索路径。因为余额都是 12 元,所以后续的求解路径和结果完全相同。

    在这个硬币求解问题中,当余额相同的时候,最优解是确定的。那么,如果能够避免相同余额下的重复搜索过程,那么算法执行速度是不是可以加快了?这就是重叠子问题

    你可以把求解 12 元的硬币数量理解成求解 25 元的硬币数量的一个子问题。在求解 25 元硬币过程中,会有很多种情况都要求解 12 元硬币的最优解。这类会出现重复求解的子问题称之为重叠子问题。

    总结与升华

    贪心算法只能解决局部最优问题,而我们的最终目标是解决整体最优问题(即最优解)。

    自然地,枚举是获得最优解的理想方法。而递归可以帮助我们获得所有可能答案的组合。递归形式的求解几乎就是简单地把题设中的函数表达式照搬过来,它相较于迭代来说更直观,且易于理解。

    但暴力递归有着十分明显的缺陷,它存在性能低下、可读性低和调试困难等问题。为此,我们可以进行剪枝:

    • 利用预设条件减少搜索路径,优化最优组合搜索方案(硬币的优化);
    • 利用重叠子问题,避免重叠子问题的计算。
    展开全文
  • 接着利用Shi--Tomasi算法对建筑物局部无法规整的复杂轮廓区域进行特征角点提取、匹配、排序剔除,进一步对边缘特征点进行有序连接重构,实现轮廓深度的优化,最终提高边缘表达准确度和提取精度。通过对多幅遥感影像...
  • 改进算法将具有信息定向流动的闭环拓扑结构星型拓扑结构或四边形拓扑结构相结合,促使粒子在前期寻优过程中具有较高的多样性,保证搜索的广度,而在后期满足粒子群的整体收敛性,保证寻优的精度.同时,将布谷鸟搜索算法...
  • 通过构造斗杆典型工况下的三维参数化优化模型,利用Ansys有限元方法分析对其危险部位的应力...结果表明优化模型可行,优化后的断面结构可以在保证斗杆重量最小化的同时,改善局部的应力分布状况,提高斗杆整体的机械性能。
  • 算法在传统禁忌搜索算法的基础上, 利用以长期记忆为基础的多元化扰动策略和块变异算子来扩大贴片机贴装顺序优化搜索空间, 并结合局部下降搜索策略优化喂料器分配, 最终实现贴片机贴装整体优化. 仿真实验表明, 改进...
  • 数学建模之优化模型详解

    千次阅读 2022-03-09 21:51:29
    全文共8090个字,码字总结不易,老铁们来...企业推出新产品要综合考虑成本市场吸引力,对资金进行优化配置,等等。 这些问题都是“最优化问题”,也是数学建模中的典型问题,解决最优化问题的数学方法就是“最优化..

    全文共8090个字,码字总结不易,老铁们来个三连:点赞、关注、评论
    作者:[左手の明天]
     原创不易,转载请联系作者并注明出处
    版权声明:本文为博主原创文章,遵循 CC 4.0 BY-SA 版权协议,转载请附上原文出处链接和本声明。

    “优化”是生活中经常使用的词:坐出租车时希望司机不绕弯路、走优化路线;逛超市时考虑各种优惠活动,希望获得最大优惠;企业推出新产品要综合考虑成本与市场吸引力,对资金进行优化配置,等等。 这些问题都是“最优化问题”,也是数学建模中的典型问题,解决最优化问题的数学方法就是“最优化方法”。

    最优化方法的出发点是系统思维,最优化方法的基本思路是在一定的约束条件下,保证各方面资源的合理分配, 最大限度地提升系统某一性能或系统整体性能,最终实现最理想结果。运用最优化方法建立并求解数学模型,主要包括以下步骤:

    (1)明确目标,分析问题背景,确定约束条件,搜集全面的客观数据和信息;
    (2)建立数学模型,构建变量之间的数学关系,设立目标函数;
    (3)分析数学模型,综合选择最适合该模型的优化方法;
    (4)求解模型,通常借助计算机和数学分析软件完成;
    (5)对最优解进行检验和实施。

    目录

    数学规划的一般模型

    MATLAB 求解优化问题的主要函数

    模型及基本函数

    优化函数的输入变量

    优化函数的输出变量

    无约束最优化问题

    数学描述

    解析解法和图解法

    数值解法

    全局最优解和局部最优解

    带约束最优化问题

    线性规划问题

    情况一

    情况二

    二次规划问题

    非线性规划问题

    定义

    求解算法1:间接法

    求解算法2:直接法

    求解算法3:最速下降法(steepest descent method)

    Matlab求解步骤

    示例

    0-1规划问题

    钢管的订购与运输问题

    问题

    问题1的基本模型和解法

    总费用最小的优化问题

    基本模型:二次规划

    Floyd算法求解步骤 

    最优化方法在数学建模中的应用 

    梯度下降法

    惩罚函数法

    遗传算法

    蚁群算法


    数学规划的一般模型

    其中,x~决策变量;f(x)~目标函数;gi(x)≤0~约束条件


    MATLAB 求解优化问题的主要函数

    模型及基本函数

    优化函数的输入变量

    优化函数的输出变量


    无约束最优化问题

    数学描述

    解析解法和图解法

     举例:用解析和图解法求解下列方程

    <<syms t; 
    y=exp(-3*t)*sin(4*t+2)+4*exp(-0.5*t)*cos(2*t)-0.5;
    ezplot(y,[0 4])
    y1=diff(y);
    ezplot(y1,[0 4])
    t0=solve(y1)
    y2=diff(y1);
    b=subs(y2,t,t0)

    数值解法

     命令形式1:

    x=fminsearch(fun,x0)   %简单形式
    
    [x,f,flag,out]=fminsearch(fun,x0,opt,p1,p2,…) %一般形式

    功能:与fsolve()中的参数控制形式类似。

    注:若函数时多元的,要表达成向量的形式。

    命令形式2:

    x=fminunc(fun,x0)   %简单形式
    
    [x,f,flag,out]=fminunc(fun,x0,opt,p1,p2,…) %一般形式

    功能:与fsolve()中的参数控制形式类似。

    举例:

    >>f=inline('(x(1)^2-2*x(1))*exp(-x(1)^2-x(2)^2-x(1)*x(2))','x');
    x0=[0,0];
    ff=optimset;ff.Display='iter';
    x=fminsearch(f,x0,ff)
    
    >>x=fminunc(f,x0,ff)
    

    全局最优解和局部最优解

    一元函数极小

    X=fminbnd(fun,x1,x2)

    多元无约束极小

    X=fminunc(fun,x0) (牛顿法)
    X=fminsearch(fun,x0)

    举例1:(初值的影响力)设目标函数为

     试观察不同的初值得出的最小值。

    >> f=inline('exp(-2*t)*cos(10*t)+exp(-3*(t+2))*sin(2*t)','t');
    t0=1;[t1,f1]=fminsearch(f,t0)
    
    t1=0.92275390625000,f1=-0.15473299821860
    
    >> t0=0.1;[t2,f2]=fminsearch(f,t0)
    
    t2=0.29445312500000,f2=-0.54362463738706
    
    
    >> syms t; 
    y=exp(-2*t)*cos(10*t)+exp(-3*(t+2))*sin(2*t);
    ezplot(y,[0,2.5]); axis([0 2.5 -0.6 1])

    举例2:对边长为3米的正方形铁板,在四个角剪去相等的正方形以制成方形无盖水槽,问如何剪法使水槽的容积最大?

    建立模型: 

    设剪去的正方形的边长为x,,则水槽的容积为:

    建立无约束优化模型为:

    模型求解:

    先编写M文件如下:

    function f=myfun(x)
    f=-(3-2*x).^2*x;

    调用fminbnd:

    [x,fval]=fminbnd(@myfun,0,1.5)

    运算结果为:

    x = 0.5000,fyal =2.0000.

    即剪掉的正方形的边长为0.5米时水槽的容积最大,最大容积为2立方米。


    带约束最优化问题

    线性规划问题

    目标函数:

     约束条件:

    情况一

    目标函数:

    其中,C为价值向量

    约束条件:

     其中,b为资源向量;X为决策变量向量

    其中:

     

    命令形式1:

    [X,lag,how]=lp(C,A,b,v1,v2,x0)

    功能:

    • C,A,b的意义如矩阵表示里参数;
    • v1,v2表示决策变量的上界和下界(其维数可以小于X,但表示前几个分量的上下界);
    • x0表示初始值;X时输出最优解;
    • lag是lagrange乘子,维数等于约束条件的个数,非零的向量是起作用的约束条件;
    • how给出错误信息:infeasible(无可行解),unbounded(无界解),ok(求解成功).

     举例:

    >> c=[13,-1,5];
    A=[-1,-1,0;0,1,1];
    b=[-7,10];
    v0=[2,0,0];
    [X,lag,how]=lp(c,A,b,v0)
    

    情况二

    目标函数:

     约束条件:

    命令形式2:  

    [X,f,flag,c]=linprog(C,A,b,Aeq,Beq,xm,xM,x0,opt)

    功能:各个参数的解释如前,若各个约束条件不存在,则用空矩阵来代替。

    • x: 解
    • f: 最优值
    • flag:大于零表示求解成功,否则求解出问题
    • c:求解信息
    • x0:搜索点的初值
    • opt:最优化控制项

    举例1:

    >> c=[-2,-1,-4,-3,-1];
     A=[0 2 1 4 2;3 4 5 -1 -1];
     b=[54,62];
     Ae=[];Be=[];
     xm=[0,0,3.32,0.678,2.57];
     ff=optimset;
     ff.LargeScale='off';
     ff.TolX=1e-15;
     ff.Display='iter';
     [X,f,flag,c]=linprog(c,A,b,Ae,Be,xm,[],[],ff)
    

    举例2:某车间生产A和B两种产品,为了生产A和B,所需的原料分别为2个和3个单位,所需的工时分别为4个和2个单位。现在可以应用的原料为100个单位,工时为120个单位。每生产一台A和B分别可获得利润6元和4元。应当生产A和B各多少台能获得最大利润?

    分析:

     模型建立:

    设生产A产品x1 台,生产B产品 x2台

     模型求解:

    f=[-6,-4]';
    A=[2 3;4 2];
    B=[100;120];
    Ae=[];
    Be=[];
    xm=[0,0];
    ff=optimset;
    ff.LargeScale='off'; % 不用大规模问题求解
    ff.TolX=1e-15;
    ff.TolFun=1e-20; 
    ff.TolCon=1e-20;
    ff.Display='iter';
    [x,f_opt,key,c]=linprog(f,A,B,Ae,Be,xm,[],[],ff)

    举例3:(任务分配问题)某车间有甲、乙两台机床,可用于加工三种工件。假定这两台车床的可用台时数分别为 800 和 900,三种工件的数量分别为 400、600 和 500,且已知用三种不同车床加工单位数量不同工件所需的台时数和加工费用如下表。问怎样分配车床的加工任务,才能既满足加工工件的要求,又使加工费用最低?

    模型建立:

    设在甲车床上加工工件 1、2、3 的数量分别为 x1、x2、x3,在乙车床上加工工件 1、2、3 的数量分别为 x4、x5、x6。可建立以下线性规划模型:

    模型求解:

    f = [13 9 10 11 12 8];
    A = [0.4 1.1 1 0 0 0
    0 0 0 0.5 1.2 1.3];b = [800; 900];
    Aeq=[1 0 0 1 0 0
    0 1 0 0 1 0
    0 0 1 0 0 1];
    beq=[400 600 500];
    vlb = zeros(6,1);
    vub=[];
    [x,fval] = linprog(f,A,b,Aeq,beq,vlb,vub)

    举例4:某厂每日 8 小时的产量不低于 1800 件。为了进行质量控制,计划聘请两种不同水平的检验员。一级检验员的标准为:速度 25 件/小时,正确率 98%,计时工资 4 元/小时;二级检验员的标准为:速度 15 小时/件,正确率 95%,计时工资 3 元/小时。检验员每错检一次,工厂要损失 2 元。为使总检验费用最省,该工厂应聘一级、二级检验员各几名?

    模型建立:

    设需要一级和二级检验员的人数分别为 x1、x2 人,则应付检验员的工资为:

     因检验员错检而造成的损失为:

    故目标函数为:

    约束条件为:

     

     线性规划模型:

     

     模型求解:

    c = [40;36];
    A=[-5 -3];
    b=[-45];
    Aeq=[];
    beq=[];
    vlb = zeros(2,1);
    vub=[9;15];
    %调用 linprog 函数:
    [x,fval] = linprog(c,A,b,Aeq,beq,vlb,vub)

    结果:

    x =
    9.0000
    0.0000
    fval =360

    即只需聘用 9 个一级检验员。

    二次规划问题

    目标函数:

    约束条件: 

    命令形式:

     [X,f,flag,c]=quadprog(H,C,A,b,Aeq,Beq,xm,xM,x0,opt)

    功能:

    各个参数的解释如前,若各个约束条件不存在,则用空矩阵来代替。

    举例:

     

    >> c=[-2,-1,-4,-3,-1];
      A=[0 2 1 4 2;3 4 5 -1 -1];
      b=[54,62];
      Ae=[];Be=[];
      xm=[0,0,3.32,0.678,2.57];
      ff=optimset;
      ff.LargeScale='off';
      ff.TolX=1e-15;
      ff.Display='iter';
      [X,f,flag,c]=linprog(c,A,b,Ae,Be,xm,[],[],ff)
    

    非线性规划问题

    定义

    如果目标函数或约束条件中至少有一个是非线性函数时的最优化问题就叫做非线性规划问题.

    一般形式:

     其中

     

     是定义在 En 上的实值函数,简记:

     

    其它情况:

    求目标函数的最大值或约束条件为小于等于零的情况,都可通过取其相反数化为上述一般形式.

     

    其中X为n维变元向量,G(X)与Ceq(X)均为非线性函数组成的向量,其它变量的含义与线性规划、二次规划中相同.

    求解算法1:间接法

    在非线性最优化问题当中,如果目标函数能以解析函数表示,可行域由不等式约束确定,则可以利用目标函数和可行域的已知性质,在理论上推导出目标函数为最优值的必要条件,这种方法就称为间接法(也称为解析法) 。 一般要用到目标函数的导数。

    求解算法2:直接法

    直接法是一种数值方法。这种方法的基本思想是迭代,通过迭代产生一个点序列{ X(k) },使之逐步接近最优点。 只用到目标函数。 如黄金分割法、Fibonacci、随机搜索法。

    迭代法一般步骤

     注意:数值求解最优化问题的计算效率取决于确定搜索方向P (k)和步长的效率。

    求解算法3:最速下降法(steepest descent method)

    由法国数学家Cauchy于1847年首先提出。在每次迭代中,沿最速下降方向(负梯度方向)进行搜索,每步沿负梯度方向取最优步长,因此这种方法称为最优梯度法。

    特点:

    方法简单,只以一阶梯度的信息确定下一步的搜索方向,收敛速度慢; 越是接近极值点,收敛越慢; 它是其它许多无约束、有约束最优化方法的基础。 该法一般用于最优化开始的几步搜索。

    最速下降法算法:

    Matlab求解步骤

    用Matlab求解上述问题,基本步骤分三步:

    1. 首先建立M文件fun.m,定义目标函数F(X):

    function f=fun(X); 
    f=F(X);

     3. 建立主程序.非线性规划求解的函数是fmincon,命令的基本格式如下:       

    (1) x=fmincon(‘fun’,X0,A,b)    
    
    (2) x=fmincon(‘fun’,X0,A,b,Aeq,beq)    
    
    (3) x=fmincon(‘fun’,X0,A,b, Aeq,beq,VLB,VUB)          
    
    (4) x=fmincon(‘fun’,X0,A,b,Aeq,beq,VLB,VUB,’nonlcon’)
    
    (5) x=fmincon(‘fun’,X0,A,b,Aeq,beq,VLB,VUB,’nonlcon’,options)
    
    (6) [x,fval]= fmincon(...)
    
    (7) [x,fval,exitflag]= fmincon(...)  
    
    (8) [x,fval,exitflag,output]= fmincon(...)

    输入参数的几点说明:

    模型中如果没有A,b,Aeq,beq,VLB,VUB的限制,则以空矩阵[ ]作为参数传入;

    nonlcon:如果包含非线性等式或不等式约束,则将这些函数编写一个Matlab函数

    nonlcon就是定义这些函数的程序文件名;

    不等式约束  G(x)<=0

    等式约束     Ceq(x)=0.

    如果nonlcon=‘mycon’ ; 则myfun.m定义如下

    function [G,Ceq] = mycon(x)
    
    G= ... % 计算非线性不等式约束在点x处的函数值
    
    Ceq= ... %计算机非线性等式约束在点x处的函数值 

    示例

    例子1:

    2个不等式约束, 2个等式约束 3个决策变量x1,x2,x3 如果nonlcon以‘test’作为参数值,则程序test.m如下

    function [G,Ceq]=test(x)
    G(1)=x(1)*x(1)+x(2)*x(2)+x(3)*x(3)-100
    G(2)=60-x(1)*x(1)+10*x(3)*x(3)
    Ceq(1)=x(1)+x(2)*x(2)+x(3)- 80
    Ceq(2)=x(1)^3+x(2)*x(2)+x(3)- 80
    

    注意:

    [1] fmincon函数提供了大型优化算法和中型优化算法。默认时,若在fun函数中提供了梯度(options参数的GradObj设置为’on’),并且只有上下界存在或只有等式约束,fmincon函数将选择大型算法。当既有等式约束又有梯度约束时,使用中型算法。

    [2] fmincon函数可能会给出局部最优解,这与初值X0的选取有关。

    例子2:

     1.先建立M文件 fun.m,定义目标函数:

    function f=fun(x);
    
    f=exp(x(1)) *(4*x(1)^2+2*x(2)^2+4*x(1)*x(2)+2*x(2)+1);

     2.再建立M文件mycon.m定义非线性约束:

    function [g,ceq]=mycon(x)
    
    g=[1.5+x(1)*x(2)-x(1)-x(2);-x(1)*x(2)-10];
    
    ceq=[];

    3.主程序为:

    x0=[-1;1];
    A=[];b=[];
    Aeq=[1 1];
    beq=0; 
    vlb=[];
    vub=[];
    [x,fval]=fmincon('fun4',x0,A,b,Aeq,beq,vlb,vub,'mycon')

    4.运算结果为:

    x = -1.2250    1.2250        
    
    fval = 1.8951

    0-1规划问题

    数学描述:自变量的取值只能为0或1

    matlab解:

    X=bintprog(f,A,B,Aeq,Beq)

    小规模问题可以穷举

    举例:求解下面的0-1线性规划问题

     

    f=[-3,2,-5]; A=[1 2 -1; 1 4 1; 1 1 0; 0 4 1]; 
    B=[2;4;5;6];
    x=bintprog(f,A,B,[],[])'

    钢管的订购与运输问题

    要铺设一条A1→A2 →……→ A15的输送天然气的主管道,如图一所示(见下页)。经筛选后可以生产这种主管道钢管的钢厂有 。图中粗线表示铁路,单细线表示公路,双细线表示要铺设的管道(假设沿管道或者原来有公路,或者建有施工公路),圆圈表示火车站,每段铁路、公路和管道旁的阿拉伯数字表示里程(单位km)。

    为方便计,1km主管道钢管称为1单位钢管。一个钢厂如果承担制造这种钢管,至少需要生产
    500个单位。钢厂Si在指定期限内能生产该钢管的最大数量为Si个单位,钢管出厂销价1单位钢管为pi万元,如下表:

     1单位钢管的铁路运价如下表:

     

    1000km以上每增加1至100km运价增加5万元。公路运输费用为1单位钢管每公里0.1万元(不足
    整公里部分按整公里计算)。钢管可由铁路、公路运往铺设地点(不只是运到点A1,A2,……,A15 ,而是管道全线)。

    问题

    (1)制定钢管的订购和运输计划,使总费用最小.
    (2)分析对购运计划和总费用影响:哪个钢厂钢管销价的变化影响最大;哪个钢厂钢管产量上限的变化影响最大?
    (3)讨论管道为树形图的情形

    问题1的基本模型和解法

    总费用最小的优化问题

     

    基本模型:二次规划

    Floyd算法求解步骤 

     Floyd算法过程描述如下:

    • 1、 首先S以边集M初始化,得到所有的直接连通代价;
    • 2、 依次考虑第k个结点,对于S中的每一个S[i][j],判断是否满足:S[i][j]>S[i][k]+S[k][j],如果满足则用S[i][k]+S[k][j]代替S[i][j],此为第k步;
    • 3、 k循环取遍所有结点,算法结束时,S为最终解。

     


    最优化方法在数学建模中的应用 

    梯度下降法

    梯度下降法是经典的最优化方法之一[4],其核心思想是高等数学中的导数理论。 梯度下降法实现最优化的原理是,每次迭代更新目标函数时,都以该变量导数(即梯度)的反方向作为更新参数的方向,最终解一定会收敛于最优解。 这个原理类似于走下坡路时,总是沿着最陡峭的方向向下走,最后就一定会走到坡底。梯度下降法的实现简单, 但是求解计算时间长,因此基于梯度下降法发展了很多改进算法,包括随机梯度下降法、小批量梯度下降法等,能够有效改善计算成本高的问题。

    惩罚函数法

    惩罚函数法,指的是引入惩罚因子和惩罚函数的最优化方法[5]。 具体来说,惩罚函数的思想是:将最优化问题中的约束条件视为围墙,而迭代更新的解视为在围墙内运动的粒子,一旦粒子靠近围墙,对应的惩罚因子数值就会增大,导致惩罚函数值增大,反之,粒子远离围墙时,惩罚函数值就减小。 建立了这种惩罚机制后,在每次迭代过程中,模型为了“避免被惩罚”,逐渐趋近于约束边界,从而找到了最优解。惩罚函数法对模型的训练虽然“简单粗暴”,但是原理直观、实现门槛低,是实际工程中备受青睐的最
    优化方法。

    遗传算法

    不同于梯度下降法和惩罚函数法,遗传算法并非依据导数理论提出的算法[6],而是一种模拟生物在自然届中进化规律的一种智能算法。 自然界的生物进化遵循适者生存和优胜劣汰,即能够适应环境变化或基因变异的个体才能够参与到进化。 遗传算法的优化原理与之类似:每一次迭代时,通过计算各个个体的适应度,从中随机地选择两个个体作为父母,繁殖后代,同时诱发子代的染色体变异,重复迭代,当出现最大适应度的子代时,即认为获得了最优解,循环结束。与梯度下降法、惩罚函数法相比,遗传算法以生物进化为原型,收敛性较好,在计算精度要求时,具有计算时间少、鲁棒性高的优势。

    蚁群算法

    与遗传算法类似,蚁群算法也是受启发于生物的一种最优化方法[7]。 生物科学家发现蚂蚁经过的路上都会有一种特殊物质,并且蚁群中的蚂蚁对该物质高度敏感,由于该物质浓度越高,代表着路途长度越短,想要走“捷径”的蚁群们都会选择浓度较高的道路行走,“捷径”经过的蚂蚁越多,特殊物质的浓度就越高,物质浓度积累到一定程度,所有蚂蚁都会被吸引到最佳捷径上来,都能以最快速度找到食物了。 蚁群算法解决最优化问题,就是利用了其分布计算和信息正反馈的特点。

    展开全文
  • 本文提出了一类基于局部-整体反馈控制的多个体系统稳定性问题的框架.该问题具有以下特点:1)每个 个体有自身的动力学...2)每个个体通过优化自身的局 部目标函数其他个体耦合在一起.在此框架下,本文给出了最优的局部
  • 介绍了敏感场的计算方法,借助大型ANSYS有限元分析工具,建立电容传感器结构模型,并得到电极间敏感场的分布,为了避免粒子在全局附近振荡和粒子陷入局部最优,提出了改进PSO改进混沌搜索策略的混合新算法....
  • 实验结果表明:所提算法去雾图像的主观视觉效果较好,且对比度和色偏程度两方面客观评价指标整体优于其他对比算法。该算法能够有效去除夜间图像雾气,提高图像的对比度,恢复更多的细节信息,且颜色自然,视觉效果...
  • 将频谱分配的二进制编码转换为量子序列编码,提出一种基于量子果蝇优化的...仿真结果表明,改进的QFOA收敛速度快且跳出局部最优能力强,应用到认知无线网络频谱分配中,增加了网络资源利用率,提高了网络的整体性能。
  • 提出一种基于多代理模型的优化方法,求解混合整数规划问题....再次,通过多代理模型对候选解进行预筛选,实现粒子群算法的协同优化;最后,通过14个测试问题和一个基于数据驱动的模型参数选取问题,验证所提出方法的有效性.
  • 首先教你SQL整体优化、快速优化实施、如何读懂执行计划、如何左右执行计划这四大必杀招。整这些干嘛呢?答案是,传授一个先整体后局部的宏观解决思路,走进“道”的世界。 接下来带领大家飞翔在“术”的天空。教你...
  • 提出了一种基于压缩感知的邻域优化算法,运用压缩感知技术对高维空间目标点近邻进行压缩采样,构建“收—放”模型,自适应得到最优子空间,同时优化邻域组成元素,使得数据的整体降维效果更加稳定。通过手工流形和...
  • 【最优化】最优化的相关条件

    千次阅读 2020-11-20 11:19:41
    局部极值点的充分、必要条件的陈述证明

    最优化的相关条件


    一. 局部极值点的必要条件

    既然在前一节中我们对于局部(全局)极值点进行了讨论,那么我们也想要知道,如果一个点是极值点,那么它需要满足什么条件?

    1. 局部极小点的必要条件
    在这里插入图片描述
    若点x0是函数f(x)的局部极小值点,则有:

    • x0是f(x)的驻点——(f(x)在x0处的一阶导数为0)
    • f(x)在x0处的Hessian矩阵H(x0)是半正定的

    2. 证明
    在这里插入图片描述

    【自我小结】

    • 有关极值的相关证明经常用到极值的定义和泰勒展开
    • 反证法也是经常用到的证明思路
    • 为了构造出上面所需的定义和泰勒展开式,往往会取某个点的邻域进行相应的演算。

    3. 局部极大点的必要条件

    前面虽然都是在以局部极小点为例进行探讨,但是我们可以很容易地将结论迁移类比到局部极大值上面。
    在这里插入图片描述
    根据上图我们可知,如果一个点x0是f(x)的局部极大值点,则有:

    • f(x)在x0处的一阶导数值为0
    • f(x)在x0的Hessian矩阵是负半定的

    二. 局部极值点的充分条件

    1. 局部极小值的充分条件
    在这里插入图片描述
    若函数f(x)上的某一点x0满足:

    • x0是f(x)的驻点——(f(x)在x0处的一阶导数为0)
    • f(x)在x0处的Hessian矩阵H(x0)是正定的

    则可得到,x0就是f(x)的严格的局部极小值点(<)

    2. 证明
    在这里插入图片描述
    3. 注意事项

    ①对照充分条件和必要条件,可知充分条件比必要条件更加严苛(必要条件涵盖了充分条件)
    ②充分条件的适用范围更窄,函数需要定义在开集上,且函数需要二阶可微

    在实际问题中,函数往往是定义在闭集上,而且最优点很可能就出现在端点上。

    这个时候,我们可以先利用充分条件找到该闭集对应的开集中的最优点,再和端点上的函数值进行比较,从而确定整个闭集上的最优点。

    ③用充分条件可以直接找到极值点中的所有非奇异点,奇异点则较难寻找

    在这里插入图片描述
    注意分辨“充分条件和“必要条件”的关系

    • 满足了充分条件的点一定是非奇异极值点,不满足充分条件的点不一定不是极值点
    • 所有极值点都必须满足必要条件,满足了必要条件不一定就是极值点
    • 满足了必要条件,但不满足充分条件的那些极值点就称为奇异点

    ④对于一个定义在开集X上的可微凸函数f(x)而言,▽f(x*) = 0↔点x*是f(x)的全局极小值点

    【简要证明】
    根据上一篇博文《【最优化】最优化理论的基本概念》中的《可微函数的凸性三定理》的定理一:

    我们可以任取一个x∈X,则有f(x)-f(x*)≥▽f(x*)T·(x-x*);
    又因为根据题意,▽f(x*)T = 0,所以可得f(x)-f(x*)≥0

    从而证得x*就是f(x)的全局极小值点。


    【例】求解给定函数的局部极值点
    在这里插入图片描述
    关于极值点求解的题型要注意:

    • 求解出驻点之后,一定要用二阶导进行验证,才能明确是不是极小(大)值点
    • 只能得出该点满足必要条件,则不能下结论说该点是局部极值点
    • 不同函数在某点的一阶导和二阶导处具有相同的性质,这并不意味着该点在两个函数中是相同性质的点,可以通过计算更高阶导数进行分析

    三. 局部极值点的存在条件

    不是所有的函数都存在极值点:
    在这里插入图片描述
    1. 极值点的存在条件
    在这里插入图片描述
    总的来说,我们只需要记住在一个给定的有界区间内求解,一般都是有极值点的

    2. 讨论极值点的各类条件的原因

    本节我们讨论了极值点的充分、必要、充分必要及存在性条件。

    • 帮助理解极值求解过程
    • 当条件满足时,我们可以适时地终止一个迭代算法
    • 可以帮助优化算法
    展开全文
  • 针对粒子群算法整体上容易陷入局部最优的缺陷,将鱼群算法中的视距、拥挤度引入标准粒子群算法,提出一种改进的粒子群算法,有效提高了粒子群算法的全局收敛性。通过基准函数Sphere、Griewank、Ackley和Shekel’s ...

空空如也

空空如也

1 2 3 4 5 ... 20
收藏数 78,628
精华内容 31,451
关键字:

局部优化与整体优化