如何用Matlab等Jacobian语言进行约束优化？

Question

我在Matlab中有一个进行约束优化的代码：

这一项是根据约束和情况导出函数的目标函数：

function f=funcobj(X);
    f=[(2-3*X(1)*2+X(3)*(1-X(1))+X(4)*(5+X(1)/5));
    (3-4*X(2)+3*X(3)+2*X(4));
    (X(1)+3*X(2)-X(1).^2/2-5.5);
    (5*X(1)+2*X(2)+X(1).^2/10-10)];
end

这个是Jacobian函数

function [f0,jac]=jacobian(x);
h = 1.0e-4;
n = length(x);
jac = zeros(n,n);
f0 = funcobj(x);
for zz=1:n;
    temp = x(zz);
    x(zz)= temp + h;
    f1 = funcobj(x);
    x(zz)= temp;
    jac(:,zz) = (f1 - f0)/h;
    disp((f1 - f0)/h);
end
% disp(jac)

这个是主要代码

clc;close all;clear all
%initial
X(1,:)= [1 1 1 1]';
niter=30;tol=1e-6;
for ii=1:niter-1
    disp(X(ii,:));
    [f,dp]=jacobian(X(ii,:));
    dX=inv(dp)*f;
    X(ii+1,:)=X(ii, :)'-dX;
    fprintf('Iterasion=ti Solution=$.4f \n',ii,X(ii+1))
    if abs(X(ii+1,:)-X(ii,:))<tol;
        r=X(ii+1,:);
        disp('The Solution is convergent')
        break
    end
end
x=r(1);y=r(2);lambda_1=r(3);lambda_2=r(4);
f = (2*x)+(3*y)-(x).^3-2*(y.^2);
disp('Case 1')
disp(['x=' num2str(x) ', y = ' num2str(y),',f = ' num2str(f)])
disp(['lambda_1 = ' num2str(lambda_1), ', lambda_2 = ' num2str(lambda_2)])

当我试图将它转换为Python时，我仍然混淆了X数组以及如何用Python重写jacobian。这是我的尝试：

import numpy as np
def funcobj(z):
    f = np.array([[2-3*z[0]**2 + z[2]*(1-z[0])+z[3]*(5+z[0]/5)], [(-4*z[1]+3*z[2]+2*z[3])], [z[0]+3*z[1]-(z[2]**2)/2-5.5], [5*z[0]+2*z[1]+(z[0]**2)/10-10]])
    print(f)
    return f

def jacobian(X):
    h = 1.0e-4
    n = len(X)
    print(n)
    jac = np.zeros([4,4])
    f0 = funcobj(X)
    for i in range(0,n):
        temp = X[i]
        X[i] = temp + h
        f1 = funcobj(X)
        X[i] = temp
        #print((f1-f0)/h)
        jac[0,i] = (f1-f0)/h
    return (f0, jac)

X=np.array([[1],[1],[1],[1]])
niter=30
tol=1e-6
for i in range(0,niter):
 jacobian(X[:,i])
 if abs(X[:,i]-X[:,i-1])<tol:
     r=X[:,i]
     print('The Solution is convergent')
     break

如何修正这段代码？我仍然知道Python中的错误

Answer 1

您的funcobj返回一个形状为(n,1)而不是(n,)的np.ndarray。注意，与matlab相反，在numpy中，前者对应于矩阵，后者对应于向量。接下来，在行jac[0, i] = (f1-f0)/h中，您尝试将一个np.ndarray分配给单个矩阵元素。它应该是jac[:, i]。还请注意，range在默认情况下从0开始，因为python有基于0的索引。

代码：

def funcobj(z):
    f1 = 2-3*z[0]**2 + z[2]*(1-z[0])+z[3]*(5+z[0]/5)
    f2 = (-4*z[1]+3*z[2]+2*z[3])
    f3 = z[0]+3*z[1]-(z[2]**2)/2-5.5
    f4 = 5*z[0]+2*z[1]+(z[0]**2)/10-10
    return np.array((f1, f2, f3, f4))

def jacobian(X):
    h = 1.0e-4
    n = len(X)
    jac = np.zeros([4,4])
    f0 = funcobj(X)
    for i in range(n):
        temp = X[i]
        X[i] = temp + h
        f1 = funcobj(X)
        X[i] = temp
        #print((f1-f0)/h)
        jac[:,i] = (f1-f0)/h
    return (f0, jac)

x0 = np.ones(4)
# works as expected
print(jacobian(x0))

现在轮到您从这里开始，在python中实现主要的算法。