使用fmincon（）解决矩阵形式的特殊非线性程序

Question

我的问题是以下优化问题：

min J=X'*E*X+U'*E*U
s.t. X'-X0'-X'*D1*Q*P+X'*D1*Q*Z=0

其中X和U是2*r+1 by 1列矩阵，X0是2*r+1 by 1已知列矩阵，而E, D1, D2, P and Z是已知2*r+1 by 2*r+1矩阵。此外，其中Phi*Phi'*D2*U=Q*Phi是Phi已知列矩阵的2*r+1 by 1。

已知矩阵：X0, E, D1, D2, Z, P, Phi

    r=2;

    % X0
      X0(1)=1;
    for i=2:2*r+1
        X0(i)=0;
    end
    X0=X0';
    % P
    P1=[1/2];
    P2=zeros(1,r);
    for i=1:r
        P3(i)=(-1)/(i*pi);
    end
    P4=zeros(r,1);
    P5=zeros(r,r);
    for i=1:r
        V1(i)=1/(2*i*pi);
    end
    P6=diag(V1);
    for i=1:r
        W(i)=1/(2*i*pi);
    end
    P7=W';
    for i=1:r
        V2(i)=(-1)/(2*i*pi);
    end
    P8=diag(V2);
    P9=zeros(r,r);
    P=3*[P1 P2 P3 ; P4 P5 P6 ; P7 P8 P9];
    % D1
    M1=[1];
    M2=zeros(1,r);
    M3=zeros(1,r);
    M4=zeros(r,1);
    for i=1:r
        V4(i)=cos((2*i*pi)/3);
    end
    M5=diag(V4);
    for i=1:r
        V5(i)=sin((2*i*pi)/3);
    end
    M6=diag(V5);
    M7=zeros(r,1);
    for i=1:r
        V6(i)=-sin((2*i*pi)/3);
    end
    M8=diag(V6);
    for i=1:r
        V7(i)=cos((2*i*pi)/3);
    end
    M9=diag(V7);
    D1=[M1 M2 M3 ; M4 M5 M6 ; M7 M8 M9];

    % D2
    N1=[1];
    N2=zeros(1,r);
    N3=zeros(1,r);
    N4=zeros(r,1);
    for i=1:r
        VV4(i)=cos((2*i*pi*2)/3);
    end
    N5=diag(VV4);
    for i=1:r
        VV5(i)=sin((2*i*pi*2)/3);
    end
    N6=diag(VV5);
    N7=zeros(r,1);
    for i=1:r
        VV6(i)=-sin((2*i*pi*2)/3);
    end
    N8=diag(VV6);
    for i=1:r
        VV7(i)=cos((2*i*pi*2)/3);
    end
    N9=diag(VV7);
    D2=[N1 N2 N3 ; N4 N5 N6 ; N7 N8 N9];

    % Z
    Z1=[1];
    Z2=zeros(1,2*r);
    for i=1:r
        Z3(i)=(3/(2*i*pi))*sin((2*i*pi)/3);
    end
    Z3=Z3';
    Z4=zeros(r,2*r);
    for i=1:r
        Z5(i)=(3/(2*i*pi))*(1-cos((2*i*pi)/3));
    end
    Z5=Z5';
    Z6=zeros(r,2*r);
    Z=[Z1 Z2 ; Z3 Z4 ;Z5 Z6];

     % E
    V3(1)=2;
    for i=2:2*r+1
        V3(i)=1;
    end
    E=diag(V3);

    % PHi
     PH1=@(x) arrayfun(@(i)cos(i*pi*x),1:r);
     PH2=@(x) arrayfun(@(i)sin(i*pi*x),1:r);
     PH = @(X) [1, PH1(X), PH2(X)];

谢谢。

Answer 1

我的尝试：

对象功能文件名：objective_function.m

function objective = objective_function(input,r, E)

X = input(1:2*r+1);
U = input(2*r+2:end);
objective = (X')*E*X +(U')*E*U;
end

约束函数文件名：constraint.m

function [inequality, equality1, equality2] = constraint(input, r,  X0, Z, Phi, Q,  P, D1, D2)
r=3;
X = input(1:2*r+1);
U = input(2*r+2:end);
% No inequality constraint 
inequality = [];
equality1 = X'- X0'-(X')*D1*Q*P+(X')*D*Q*Z;
equality2 = Phi*Phi'*D2*U=Q*Phi;
end

优化文件名：main.m

    clear all
    clc
    r=2;

    % X0
      X0(1)=1;
    for i=2:2*r+1
        X0(i)=0;
    end
    X0=X0';
    % P
    P1=[1/2];
    P2=zeros(1,r);
    for i=1:r
        P3(i)=(-1)/(i*pi);
    end
    P4=zeros(r,1);
    P5=zeros(r,r);
    for i=1:r
        V1(i)=1/(2*i*pi);
    end
    P6=diag(V1);
    for i=1:r
        W(i)=1/(2*i*pi);
    end
    P7=W';
    for i=1:r
        V2(i)=(-1)/(2*i*pi);
    end
    P8=diag(V2);
    P9=zeros(r,r);
    P=3*[P1 P2 P3 ; P4 P5 P6 ; P7 P8 P9];
    % D1
    M1=[1];
    M2=zeros(1,r);
    M3=zeros(1,r);
    M4=zeros(r,1);
    for i=1:r
        V4(i)=cos((2*i*pi)/3);
    end
    M5=diag(V4);
    for i=1:r
        V5(i)=sin((2*i*pi)/3);
    end
    M6=diag(V5);
    M7=zeros(r,1);
    for i=1:r
        V6(i)=-sin((2*i*pi)/3);
    end
    M8=diag(V6);
    for i=1:r
        V7(i)=cos((2*i*pi)/3);
    end
    M9=diag(V7);
    D1=[M1 M2 M3 ; M4 M5 M6 ; M7 M8 M9];

    % D2
    N1=[1];
    N2=zeros(1,r);
    N3=zeros(1,r);
    N4=zeros(r,1);
    for i=1:r
        VV4(i)=cos((2*i*pi*2)/3);
    end
    N5=diag(VV4);
    for i=1:r
        VV5(i)=sin((2*i*pi*2)/3);
    end
    N6=diag(VV5);
    N7=zeros(r,1);
    for i=1:r
        VV6(i)=-sin((2*i*pi*2)/3);
    end
    N8=diag(VV6);
    for i=1:r
        VV7(i)=cos((2*i*pi*2)/3);
    end
    N9=diag(VV7);
    D2=[N1 N2 N3 ; N4 N5 N6 ; N7 N8 N9];

    % Z
    Z1=[1];
    Z2=zeros(1,2*r);
    for i=1:r
        Z3(i)=(3/(2*i*pi))*sin((2*i*pi)/3);
    end
    Z3=Z3';
    Z4=zeros(r,2*r);
    for i=1:r
        Z5(i)=(3/(2*i*pi))*(1-cos((2*i*pi)/3));
    end
    Z5=Z5';
    Z6=zeros(r,2*r);
    Z=[Z1 Z2 ; Z3 Z4 ;Z5 Z6];

     % E
    V3(1)=2;
    for i=2:2*r+1
        V3(i)=1;
    end
    E=diag(V3);

    % PHi
     PH1=@(x) arrayfun(@(i)cos(i*pi*x),1:r);
     PH2=@(x) arrayfun(@(i)sin(i*pi*x),1:r);
     PH = @(X) [1, PH1(X), PH2(X)];

     A=[];
    b=[];
    Aeq=[];
    beq=[];
    lb=[];
    ub=[];

    initial=ones(4*r+2,1);

% Objective function
J =@(decision_variable)objective_function(decision_variable, r, E);

% Constraint
equality1 = @(decision_variable)constraint(decision_variable, r,  X0, Z, Phi, Q,  P, D1, D2);
equality2 = @(decision_variable)constraint(decision_variable, r,  X0, Z, Phi, Q,  P, D1, D2);


solution = fmincon(J,initial,A,b,Aeq,beq,lb,ub,equality1, equality2);

% X and U extraction 
 X_sol = solution(1:2*r+1);
 U_sol = solution(2*r + 2:end);