开发者

解读python cvxpy下SDP问题编程

开发者 https://www.devze.com 2022-12-18 09:32 出处:网络 作者: zy123457
目录python cvxpy下SDP问题编程算法如下附上代码总结python cvxpy下SDP问题编程 最近在做定位算法的复现问题,遇到了
目录
  • python cvxpy下SDP问题编程
    • 算法如下
    • 附上代码
  • 总结

    python cvxpy下SDP问题编程

    最近在做定位算法的复现问题,遇到了

    Source Localization in Wireless Sensor Networks From Signal Time-of-Arrival Measurements

    里面的一个半正定优化算法,因此选用cvxppythony库实现。

    官方文档[cvxpy]的例程复现的算法都很简单,因此对该问题的借鉴意义不大。

    算法如下

    解读python cvxpy下SDP问题编程

    对我而言,首先的难度就是拼接矩阵后的半正定约束条件,起初是另设立两个矩阵变量,然后按部就班的增加限制条件。但最后求得的数据千奇百怪,与预测位置没有任何关系。

    后来不js断尝试更改约束限制的表达形式,但均无效果。

    后来输出了每个变量的值查看,发现Q元素的物理意义为预测距离的平方,但是求出来的Q矩阵元素往往极大,因此擅自添加了一个约束条件,限制Q的最大元素在预测距离平方的量级上,完美解决问题。

    附上代码

    class Program_t:
    
        def __init__(self,bt):
            self.BT = bt
            self.BT_x = [b[0] for b in self.BT]
            self.BT_y 编程客栈= [b[1] for b in self.BT]
            self.T=[b[2] for b in self.BT]
            self.number = len(bt)
            
        def LS_steps(self):
            num = len(self.T)
            up_control = 2*max(self.T)**2#限制最大元素量级
            Q = cp.Variable((num,num))#待求变量
            tao = cp.Variable((num,1))#生成矩阵形式后面才可以拼接
            y_ = cp.Variable((2,1))
            y_s = cp.Variable((1,1))#矩阵形式用于拼接
            yita = 0.000005*sum(self.T) / num#论文编程客栈给出的参数选择,可更改常数
            G = np.eye(num)-np.ones((num,num))
            t = np.array([self.T]).T
            expr1 = cp.trace((cp.transpose(G)) @ G @ (Q- cp.multiply(2,t @ (cp.transpose(tao)))+t @ (cp.transpose(t))))
            expr2 = yita*cp.sum(Q)
            exppythonr = expr1+ exp开发者_Python开发r2#目标函数
            Q_ = cp.bmat([[Q,tao],[cp.transpose(tao),[[1]]]])#拼接矩阵
            Y = cp.bmat([[np.eye(2),y_],[cp.transpose(y_),y_s]])#拼接矩阵
            constraints = [Q_ >> 0, Y >> 0, cp.max(Q)<=up_control]#限制条件Q半正定,Y半正定,Q最大元素小于上限(这个约束非常重要,是我自己加上去的)
            for i in range(num):
                X = np.array([self.BT_x[i],self.BT_y[i],-1]).T
                constraints += [Q[i, i] == cp.transpose(X) @ Y @ X]#约束条件
                for j in range(i+1,num):
                    X_j = np.array([self.BT_x[j], self.BT_y[j],-1]).T
                    constraints += [Q[i, j] >= cp.abs(cp.transpose(X) @ Y @ X_j)]#约束条件
            obj = cp.Minimize(expr)
            prob = cp.Problem(obj, constraints)
            prob.solve()
            position = y_.value
            print(expr1.value)#输出值
            print(expr2.value)
            print(prob.value)#输出值
            print(prob.status)#输出状态
            print(position)
            return position
           	
    

    总结

    1.理论算法与编程实现永远不等,不能轻易照搬,具体实现过程中要结合实际情况进行考虑,当求得的结果与预计相差很多时,可以尝试增加数值约束,因为计算机仿真只是近似,不是理论上的完美条件。

    2.编程实现调用库时,最好按照库的标准写,如本例中矩阵点乘可以用numpy 的dot或者cvxpy的@,以及转置的.T和cp.transpose().但是dot有时会产生意想不到的情况,平白增加工作量。

    3.复现算法时必须要对算法有深入理解,否则难以发现问题所在。

    4.不要轻易怀疑工具包的问题,经过大量使用的工具包一定比你的感觉可靠。

    以上为个人经验,希望能给大家一个参考,也希望大家多多支持我们。

    0

    精彩评论

    暂无评论...
    验证码 换一张
    取 消

    关注公众号