解讀python?cvxpy下SDP問題編程

更新時間：2022年12月17日 12:04:26 作者：zy123457

這篇文章主要介紹了解讀python?cvxpy下SDP問題編程，具有很好的參考價值，希望對大家有所幫助。如有錯誤或未考慮完全的地方，望不吝賜教

python cvxpy下SDP問題編程

最近在做定位算法的復現(xiàn)問題，遇到了

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

里面的一個半正定優(yōu)化算法，因此選用cvxpy庫實現(xiàn)。

官方文檔[cvxpy]的例程復現(xiàn)的算法都很簡單，因此對該問題的借鑒意義不大。

算法如下

對我而言，首先的難度就是拼接矩陣后的半正定約束條件，起初是另設立兩個矩陣變量，然后按部就班的增加限制條件。但最后求得的數(shù)據(jù)千奇百怪，與預測位置沒有任何關系。

后來不斷嘗試更改約束限制的表達形式，但均無效果。

后來輸出了每個變量的值查看，發(fā)現(xiàn)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#論文給出的參數(shù)選擇，可更改常數(shù)
        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)
        expr = expr1+ expr2#目標函數(shù)
        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)#輸出狀態(tài)
        print(position)
        return position