How to perform Constrained Optimization by Jacobian like Matlab?

user_18541354 0
Viewed 60

I have a code for doing constrained optimization in Matlab:

this one is for objective function from deriving a function according to the constraint and the case:

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

and this one is the Jacobian function

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)

and this one is the main code

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)])

When I try to convert it to Python, I still confused with the X array and how to rewrite jacobian in Python. This is my attempt:

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

How to fix this code? I still get the error in Python

matlab numpy optimization python
1 Answers

Your funcobj returns a np.ndarray with shape (n,1) instead of (n,). Note that contrary to matlab, in numpy, the former corresponds to a matrix while the latter corresponds to a vector. Next, in the line jac[0, i] = (f1-f0)/h you are trying to assign a np.ndarray to a single matrix element. It should be jac[:, i] instead. Note also that range starts by default at 0 since python has 0-based indexing.

In code:

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))

Now it's your turn to go on from here and implement the main algorithm in python.

user_47455290
Related