目录
1.简介
通过采用共轭梯度类似的解题步骤,用python编程实现,做到考虑问题minfx=12xTAx+bTx+c(其中 A为对称矩阵,c为常数),输入初始区间A,b,c,以及初始点,输出最优解。
2.程序设计思路
共轭梯度法的步骤:
考虑问题minfx=12xTAx+bTx+c(其中 A为对称矩阵,c为常数)
- 任意给定一个初始点x(1),置k=1
- 计算出目标函数fx在这点的梯度gk=∇fxk,若gk=0,则停止计算,得x=xk;否则,进行下一步
- 令dk=-gk+βk-1dk-1,当k=1时,βk-1=0,当k>1时,βk=gk+12gk2
- 计算λk=-g1Td(1)d1TAd1,x(k+1)=x(k)+λkd(k)
- 若k=n,则停止计算,得x=xk+1;否则,置k=k+1,返回步骤(2)
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
else:
d[p]=-g[p]+np.dot(h[0],d[p-1])
print("d=="+str(d[p]))
def caculateβ(p):
global h
#h=np.dot(np.dot(d[p],matrix),g[p+1])/np.dot(np.dot(d[p],matrix),d[p])
h=np.dot(g[p],g[p])/np.dot(g[p-1],g[p-1])
h = np.array(h)
h = h.reshape(-1)
print("β="+str(h[0]))
def caculateλ(p):
global e,d
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)
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
x=[]
getinput()
g=np.zeros((5000,n))
for i in range(2):
calulateg(i)
i = 0
j = 0
while judge(i):
print("第" + str(i + 1) + "次迭代")
calulateg(i)
if i!=0:
caculateβ(i)
caculated(i)
caculateλ(i)
i+=1
print("最优解:")
for l in range(m):
print(f'x{l+1}='+str(format(e[l], '.2f')))
if __name__ == '__main__':
main()
4.测试样例
四元:
五元: