目录
1.简介
通过采用拟牛顿法类似的解题步骤,用python编程实现,做到考虑问题minfx=12xTAx+bTx+c(其中 A为对称矩阵,c为常数),输入初始区间A,b,c,以及初始点,输出最优解。
2.程序编写思路
考虑问题minfx=12xTAx+bTx+c(其中 A为对称矩阵,c为常数)
- 任意给定一个初始点x(1)
- 置H1=In(单位矩阵),计算出目标函数fx在x1的梯度gk=∇fxk,置k=1
- 令dk=-Hkgk-1,Hk+1=Hk+p(k)p(k)Tp(k)Tq(k)-Hkq(k)q(k)THkq(k)THkq(k),p(k)=x(k+1)-xk=λkdk, q(k)=g(k+1)-gk
- 计算λk=-g1Td(1)d1TAd1,x(k+1)=x(k)+λkd(k)
- 若gk=0,则停止计算,得x=xk;否则,置k=k+1,返回(3)
3.程序代码
'''
拟牛顿法
'''
import numpy as np
from fractions import Fraction as f
def getinput():
global matrix,b,c,m,n,e
string=input('''求解的函数 min 1/2XT*A*X+bXT+C
请输入对称正定矩阵A,以;作为每一行的分隔:
''')
a = [list(map(eval,row.split())) for row in string.split(';')]
matrix = np.mat(a)
m, n = matrix.shape
string2 = input('''请输入B:''')
b=list(map(eval,string2.split(' ')))
string3=input('''请输入C:''')
c=int(string3)
string4=input('''取初始点,例如:(5,4),则输入“5 4”''')
e=list(map(eval,string4.split(' ')))
e=np.array(e)
x=[]
print('\n输入的待求解函数为')
for i in range(m):
for j in range(n):
if matrix[i,j]!=0:
x.append(f'{matrix[i,j]}*x_{i+1}*x_{j+1}')
for i in range(len(b)):
if b[i]!=0:
x.append(f'{b[i]}*x_{i+1}')
if c!=0:
x.append(f'{c}')
print('min ' + ' + '.join(x))
# print('\n初始点为:['+str(e[0])+','+str(e[1])+']')
def calulateg(p):
global g
for i in range(m):
for j in range(n):
g[p][i] = 0
for i in range(m):
for j in range(n):
# print(i)
g[p][i] += 2 * matrix[i, j] * e[j]
g[p][i] +=b[i]
def caculated(p):
global d
if p==0:
d = -g
d[p]=np.dot(-H,g[p])
print("g="+str(g[p]))
print("d="+str(d[p]))
def caculateH(p):
global H,s,f,q
if p==0:
s=H
else:
f=np.mat(f)
q=np.mat(q)
H = s + np.dot(f.T, f) / np.dot(f, q.T) - np.dot(np.dot(np.dot(s, q.T), q), s) / np.dot(np.dot(q, s), q.T)
s=H
def caculatepq(p):
global h,w,f,q
w=np.array(w)
f = np.dot(w,d[p])
q = g[p] - g[p-1]
# print("p=")
# print(f)
# print("q=")
# print(q)
def caculateλ(p):
global e,d,w
u=np.dot(-1,np.dot(g[p],d[p]))/np.dot(np.dot(d[p],np.dot(2,matrix)),d[p])
u=np.array(u)
u=u.reshape(-1)
e=np.float64(e)
d=np.float64(d)
w=u[0]
print("λ=" + str(u[0]))
for i in range(len(e)):
e[i]+=u[0]*d[p][i]
def judge(p):
if np.dot(g[p],g[p])==0:
return 0
else:
return 1
def main():
global g,H,v
x=[]
getinput()
g=np.zeros((5000,n))
H=np.identity(m)
for i in range(2):
calulateg(i)
i = 0
j = 0
while judge(i):
print("第" + str(i + 1) + "次迭代")
calulateg(i)
if i!=0:
caculatepq(i)
caculateH(i)
caculated(i)
caculateλ(i)
i+=1
print("最优解:")
for l in range(m):
print(f'x{l+1}='+str(format(e[l], '.3f')))
if __name__ == '__main__':
main()
4.测试样例
四元: