本篇文章是系列文章的第三篇,下文小编首先分享个人对二次规划定义的理解,然后举一个例题,最后使用 MindOpt Python 语言的 API 来建模以及求解 二次规划问题示例 中的问题以及展示求解的结果。
MindOpt Python、C、C++语言求解LP、MILP、QP问题系列
- Python: 线性规划LP问题、混合整数线性规划MILP问题、二次规划QP问题(本篇)
- C : 线性规划LP问题、混合整数线性规划MILP问题、二次规划QP问题
- C++ : 线性规划LP问题、混合整数线性规划MILP问题、二次规划QP问题
安装
用户可以点这里下载安装MindOpt优化求解器,免费的。找不到安装步骤点这里。
(官网https://opt.aliyun.com有更多信息等着您哟!)
二次规划定义
在前文线性规划问题示例中,讲述到线性规划在我个人认为是在线性的目标和约束中,找出一个最优解。而本文的二次规划,是非线性规划中的一类。具体地说,是一个线性约束的、二次规划问题,就是优化(最小化或最大化)二次函数目标的问题。
关于优化的类别,有很多,比如MindOpt的案例广场的标签里面提到的问题标签,就列出了常见的数学规划的类型。其中关于变量、约束、目标这建模三要素,进行罗列:
- 关于变量:取值有连续的,有整数的,还有更特殊的二进制(0或1)的,
- 关于约束和目标:一般用变量的函数变换来表达,其中约束再增加它函数的取值范围。
- 当函数是变量的线性关系时,比如x的1次方相加,我们称呼为线性约束、线性的目标。(如果变量也是连续的,这个就是线性规划问题啦。)
- 当函数是变量的是二次关系的时候,比如函数中有 x的2次方项。我们称呼为二次约束,或二次目标。
- 函数还会有凸函数和非凸函数,数学里面都代表不同的特性,大家可以再多去查阅材料。
本文主要讲 凸二次规划,Convex Quadratic Programming。
数学形式下的二次规划问题:
公式参考自MindOpt文档:https://solver.damo.alibaba.com/doc/html/model/qp/quadratic%20problem.html
案例
讲一个简单的例子,使用二次规划方法优化汽车轨迹,自动化驾驶车辆行驶在道路比较狭窄的路径上,还有其他障碍物阻碍的情况下,如果需要快速通过的话,我们需要暂时借用相邻车道通过,这个情况需要考虑自身车辆的情况、交通规则、保障远离障碍物距离的信息,然后找出一条通道。那么这个例子的解决办法是先考虑自身车辆的位置和周围障碍物,精确处理前一步可用车道,得到路径的边界,然后对路径边界进行优化(比如把车辆和障碍物之间的距离最大化,以允许车辆安全通过间隙)。
再比如,二次规划在机器人中应用,解决机械臂控制的问题。对于机器人来说,常常具备冗余关节,并且存在关节角度、角速度等限制条件,非常适合用二次规划方法来进行优化求解。
数学算例
接下来我们举一个简单的数学算例,和如何用MindOpt优化求解器进行求解。
二次规划问题示例:
Python+MindOpt代码实现
核心的用的几个APIs是:
model = MdoModel() model.set_int_attr(MDO_INT_ATTR.MIN_SENSE, 1) model.set_quadratic_elements([ x[0], x[1], x[1], x[2], x[3] ], [ x[0], x[0], x[1], x[2], x[3] ], [ 1.0, 0.5, 1.0, 1.0, 1.0 ] model.solve_prob() model.display_results()
下面是完整的例子,假设它命名为test_QP.py:
# 引入python包 from mindoptpy import * if __name__ == "__main__": MDO_INFINITY = MdoModel.get_infinity() # Step 1. 创建模型并更改参数。 model = MdoModel() try: # Step 2. 输入模型 # 改为最小化问题。 # 通过 mindoptpy.MdoModel.set_int_attr() 将目标函数设置为 最小化 model.set_int_attr(MDO_INT_ATTR.MIN_SENSE, 1) # 添加变量。 # 调用 mindoptpy.MdoModel.add_var() 来添加四个优化变量, # 定义其下界、上界、名称和类型 x = [] x.append(model.add_var(0.0, 10.0, 1.0, None, "x0", False)) x.append(model.add_var(0.0, MDO_INFINITY, 1.0, None, "x1", False)) x.append(model.add_var(0.0, MDO_INFINITY, 1.0, None, "x2", False)) x.append(model.add_var(0.0, MDO_INFINITY, 1.0, None, "x3", False)) # 添加约束。 # 注意这里的非零元素是按行顺序输入的。 model.add_cons(1.0, MDO_INFINITY, 1.0 * x[0] + 1.0 * x[1] + 2.0 * x[2] + 3.0 * x[3], "c0") model.add_cons(1.0, 1.0, 1.0 * x[0] - 1.0 * x[2] + 6.0 * x[3], "c1") # 添加二次目标矩阵 Q。 # # 注意。 # 1. 目标函数定义为c^Tx + 1/2 x^TQx,其中Q以坐标格式存储。 # 2. Q 将在内部缩放 1/2。 # 3. 为了保证Q的对称性,用户只需要输入下三角部分。 # # Q = [ 1.0 0.5 0 0 ] # [ 0.5 1.0 0 0 ] # [ 0.0 0.0 1.0 0 ] # [ 0 0 0 1.0 ] # 调用 mindoptpy.MdoModel.set_quadratic_elements() 来设置目标的二次项系数 model.set_quadratic_elements([ x[0], x[1], x[1], x[2], x[3] ], [ x[0], x[0], x[1], x[2], x[3] ], [ 1.0, 0.5, 1.0, 1.0, 1.0 ]) # Step 3. 解决问题并填充结果。 # 调用 mindoptpy.MdoModel.solve_prob() 求解优化问题, # 并用 mindoptpy.MdoModel.display_results() 来查看优化结果 model.solve_prob() model.display_results() # 调用 mindoptpy.MdoModel.get_status() 来检查求解器的优化状态, # 并通过 mindoptpy.MdoModel.get_real_attr() 和 # mindoptpy.MdoVar.get_real_attr() 来获取目标值和最优解。 status_code, status_msg = model.get_status() if status_msg == "OPTIMAL": print("Optimizer terminated with an OPTIMAL status (code {0}).".format(status_code)) print("Primal objective : {0}".format(round(model.get_real_attr(MDO_REAL_ATTR.PRIMAL_OBJ_VAL), 2))) for curr_x in x: print(" - x[{0}] : {1}".format(curr_x.get_index(), round(curr_x.get_real_attr(MDO_REAL_ATTR.PRIMAL_SOLN), 2))) else: print("Optimizer terminated with a(n) {0} status (code {1}).".format(status_msg, status_code)) except MdoError as e: print("Received Mindopt exception.") print(" - Code : {}".format(e.code)) print(" - Reason : {}".format(e.message)) except Exception as e: print("Received exception.") print(" - Reason : {}".format(e)) finally: # Step 4. 释放模型。 model.free_mdl()
MindOpt求解的结果
运行代码,比如在命令行中执行python test_QP.py,得到的结果如下所示。(#号后面是小编的注解)
Model summary. - Num. variables : 4 - Num. constraints : 2 - Num. nonzeros : 7 - Bound range : [1.0e+00,1.0e+01] #限制范围 - Objective range : [1.0e+00,1.0e+00] #目标范围 - Quad. obj. range : [5.0e-01,1.0e+00] #obj范围 - Matrix range : [1.0e+00,6.0e+00] #矩阵范围 Presolver started. Presolver terminated. Time : 0.000s Interior point method started. Iter PrimObj DualObj PrimFea DualFea GapFea Mu Time 0 +5.21950421e+01 -5.93593455e+01 1.3e+00 8.0e-01 2.1e+00 1.5e+01 0.01s 1 +5.75093325e+00 -3.28624247e+00 3.2e-02 2.0e-02 2.8e+00 1.5e+00 0.01s 2 +1.19681205e+00 +1.03397025e-04 8.1e-04 3.7e-03 1.2e+00 2.0e-01 0.01s 3 +6.52164783e-01 +3.52420863e-01 1.7e-04 3.7e-03 3.0e-01 4.9e-02 0.01s 4 +4.65540318e-01 +4.35143347e-01 4.2e-06 9.3e-05 3.0e-02 5.1e-03 0.01s 5 +4.40907312e-01 +4.39861230e-01 1.0e-07 2.3e-06 1.0e-03 1.7e-04 0.01s 6 +4.40022716e-01 +4.39996554e-01 2.6e-09 5.8e-08 2.6e-05 4.4e-06 0.01s 7 +4.40000569e-01 +4.39999914e-01 6.5e-11 1.5e-09 6.6e-07 1.1e-07 0.01s 8 +4.40000014e-01 +4.39999998e-01 1.6e-12 3.7e-11 1.6e-08 2.7e-09 0.01s 9 +4.40000000e-01 +4.40000000e-01 4.1e-14 9.1e-13 4.1e-10 6.9e-11 0.01s Terminated. - Method : Interior point method. #内点法 - Primal objective : 4.3999999966807E-01 - Dual objective : 4.3999999996074E-01 - Num. threads : 1 - Num. iterations : 9 - Solver details : Solver terminated with a primal/dual optimal status. Interior point method terminated. Time : 0.010s Optimizer terminated with an OPTIMAL status (code 1). Primal objective : 0.44 - x[0] : 0.0 - x[1] : 0.0 - x[2] : 0.2 # 决策变量的最佳取值 - x[3] : 0.2 Optimizer summary. - Optimizer used : Interior point method - Optimizer status : OPTIMAL - Total time : 0.010s Solution summary. Primal solution #目标函数最优解 - Objective : 4.3999999967e-01
联系我们
钉钉:damodi
邮箱地址:solver.damo@list.alibaba-inc.com