網站首頁 編程語言 正文
python cvxpy下SDP問題編程
最近在做定位算法的復現問題,遇到了
Source Localization in Wireless Sensor Networks From Signal Time-of-Arrival Measurements
里面的一個半正定優化算法,因此選用cvxpy庫實現。
官方文檔[cvxpy]的例程復現的算法都很簡單,因此對該問題的借鑒意義不大。
算法如下
對我而言,首先的難度就是拼接矩陣后的半正定約束條件,起初是另設立兩個矩陣變量,然后按部就班的增加限制條件。但最后求得的數據千奇百怪,與預測位置沒有任何關系。
后來不斷嘗試更改約束限制的表達形式,但均無效果。
后來輸出了每個變量的值查看,發現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)
expr = expr1+ expr2#目標函數
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.不要輕易懷疑工具包的問題,經過大量使用的工具包一定比你的感覺可靠。
原文鏈接:https://blog.csdn.net/zy123457/article/details/124681789
相關推薦
- 2022-04-05 老生常談Python中的Pickle庫_python
- 2022-11-10 Android中PopupWindow彈出式窗口使用方法詳解_Android
- 2022-08-02 python?GUI編程實現掃雷游戲_python
- 2022-06-16 基于Python+Matplotlib實現直方圖的繪制_python
- 2022-12-07 C++?IO設備讀寫功能實現詳解_C 語言
- 2022-07-25 C++超詳細講解內存空間分配與this指針_C 語言
- 2022-12-13 C++?Boost?Format超詳細講解_C 語言
- 2024-01-15 mybatis中@Results,@ResultMap注解使用
- 最近更新
-
- window11 系統安裝 yarn
- 超詳細win安裝深度學習環境2025年最新版(
- Linux 中運行的top命令 怎么退出?
- MySQL 中decimal 的用法? 存儲小
- get 、set 、toString 方法的使
- @Resource和 @Autowired注解
- Java基礎操作-- 運算符,流程控制 Flo
- 1. Int 和Integer 的區別,Jav
- spring @retryable不生效的一種
- Spring Security之認證信息的處理
- Spring Security之認證過濾器
- Spring Security概述快速入門
- Spring Security之配置體系
- 【SpringBoot】SpringCache
- Spring Security之基于方法配置權
- redisson分布式鎖中waittime的設
- maven:解決release錯誤:Artif
- restTemplate使用總結
- Spring Security之安全異常處理
- MybatisPlus優雅實現加密?
- Spring ioc容器與Bean的生命周期。
- 【探索SpringCloud】服務發現-Nac
- Spring Security之基于HttpR
- Redis 底層數據結構-簡單動態字符串(SD
- arthas操作spring被代理目標對象命令
- Spring中的單例模式應用詳解
- 聊聊消息隊列,發送消息的4種方式
- bootspring第三方資源配置管理
- GIT同步修改后的遠程分支