【问题标题】:Solving a special nonlinear program using fmincon()使用 fmincon() 求解一个特殊的非线性程序
【发布时间】:2019-12-02 08:41:36
【问题描述】:

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

min J=X'*E*X+U'*E*U
s.t. X'- X0'-(X')*D1*Q*P+(X')*D1*Q*Z=0,
     I*I'*(D2')*U-Q*I=0

其中XU2*r+1 by 1 列矩阵,X02*r+1 by 1 已知列矩阵和E, D1, D2, P and Z 是已知2*r+1 by 2*r+1 矩阵和I2*r+1 by 1 已知列矩阵.另外,Q 是一个矩阵,满足 I*I'*(D2')*U=Q*I.

已知矩阵:X0, E, D1, D2, Z, P, I

提前致谢。

我的尝试:

客观函数文件名:objective_function.m

function objective = objective_function(Unknown,r, E) 
Unknown = ones(2*r+3, 2*r+1);
X = Unknown(1, :);
U = Unknown(2, :);
Q = Unknown(3:end, :);
objective = (X')*E*X +(U')*E*U;
end

约束函数文件名:constraint.m

function [inequality, equality] = constraint(Unknown, r,  X0, Z, I, P, D1, D2)
Unknown = ones(2*r+3, 2*r+1);
X = Unknown(1, :);
U = Unknown(2, :);
Q = Unknown(3:end, :);
% No inequality constraint 
inequality = [];
equality = [X'- X0'-(X')*D1*Q*P+(X')*D1*Q*Z ; I*I'*(D2')*U-Q*I];
end

优化文件名:ma​​in.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
     R1=@(x) arrayfun(@(i)cos(i*pi*x),1:r);
     R2=@(x) arrayfun(@(i)sin(i*pi*x),1:r);
     R = @(t) [1, R1(t), R2(t)];
     I=R(1);
     I=I';
     

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

      initial=ones(2*r+3, 2*r+1);
   
% Objective function
J =@(decision_variable)objective_function(decision_variable, r, E);

% Constraint
equality = @(decision_variable)constraint(decision_variable, r, X0, Z, I, P, D1, D2);



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

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



 

【问题讨论】:

  • 您的问题是什么?你试图做什么,有什么问题?你的程序是否抛出了一些错误?您在优化中寻找哪些变量?什么是 Phi 以及它在您的优化问题中的什么位置?此外,如果你能提供你试图解决的实际数学问题以及代码本身,那将是有益的。
  • @Thales 我的变量是矩阵XU。实际上,我不知道如何在问题的约束部分定义Phi*Phi'*D2*U=Q*Phi,其中Phi 是上面定义的t 中的一个函数。我可以这样定义吗:function [inequality, equality1,equality2] = constraint(input, r, X0, E, D1, D2, Z, P, Phi) X = input(1:2*r+1); U = input(2*r+2:end); % No inequality constraint inequality = []; equality = X'-X0'-X'*D1*Q*P+X'*D1*Q*Z=0; equality2=Phi*Phi'*D2*U-Q*Phi; end
  • @Adam 你能帮我解决我的问题吗?
  • @Adam 我的意思是 Phi = @(t) [1, Phi1(t), Phi2(t)]。我在上面编辑过。我们可以把 Phi=Phi(1)。我已经编辑了我们需要的所有数学背景。
  • @Adam 当然。你说的是对的。我会努力做到的。谢谢。

标签: matlab nonlinear-optimization


【解决方案1】:

UX 是列向量描述如下r = 1

X = [x11; x12; x13] --> X is 3 by 1

U = [u11; u12; u13] --> Y is 3 by 1

Q = [ q11 q12 q13; 
      q21 q22 q23; 
      q31 q32 q33] ---> Q is 3 by 3

Unknown = [ x11 x12 x13; 
            u11 u12 u13; 
            q11 q12 q13; 
            q21 q22 q23; 
            q31 q32 q33]--> Unknown is 5 by 3

Unknown(1, :) is [ x11 x12 x13] --> 1 by 3 --> X = transpose(Unknown(1, :))

% Variable identification
---> X = Unknown(1, :)'
---> U = Unknown(2, :)'
---> Q = Unknown(3:end, :)

constraint.m 连接时存在维度问题

equality = [equality1; equality2]

equality1 = X'-X0'-(X')*D1*Q*P+(X')*D1*Q*Z

equality1_dimension = 1*3-1*3-(1*3)*(3*3)*(3*3)*(3*3)+(1*3)*(3*3)*(3*3)*(3*3) = 1*3

equality2 = I*I'*(D2')*U-Q*I

equality2_dimension = (3*1)*(1*3)*(3*3)*(3*1) - (3*3)*(3*1) = 3*1
  • 两个等式应该有相同的列数。
  • 只需转置其中一个即可。
equality = [equality1; equality2']

简而言之,更新constraint.mobjective_function.m,但请先阅读以了解原因。

目标函数文件名:objective_function.m

function objective = objective_function(Unknown,r, E) 
X = Unknown(1, :)';
U = Unknown(2, :)';
Q = Unknown(3:end, :);
objective = (X')*E*X +(U')*E*U;
end

约束函数文件名:constraint.m

function [inequality, equality] = constraint(Unknown, r,  X0, Z, I, P, D1, D2)
X = Unknown(1, :)';
U = Unknown(2, :)';
Q = Unknown(3:end, :);
% No inequality constraint 
inequality = [];
equality = [X'- X0'-(X')*D1*Q*P+(X')*D1*Q*Z ; (I*I'*(D2')*U-Q*I)'];
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
     R1=@(x) arrayfun(@(i)cos(i*pi*x),1:r);
     R2=@(x) arrayfun(@(i)sin(i*pi*x),1:r);
     R = @(t) [1, R1(t), R2(t)];
     I=R(1);
     I=I';


     A=[];
    b=[];
    Aeq=[];
    beq=[];
    lb=[];
    ub=[];
% intial contains X, U and Q
initial = ones(2*r+3, 2*r+1);

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

% Constraint
equality = @(decision_variable)constraint(decision_variable, r, X0, Z, I, P, D1, D2);



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

% X U and Q extraction 
 X_sol = solution(1, :)
 U_sol = solution(2, :)
 Q_sol = solution(3:end, :)

解决方案

X_sol =

    0.2905   -0.2881   -0.2905    0.3984   -0.2893


U_sol =

   -0.0097   -0.0097    0.0097    0.0167    0.0167


Q_sol =

    2.6741    6.1937    3.4713    1.9752    4.3157
    0.1062    0.7722    0.7143    0.4598    0.6802
    1.7584    3.6335    1.8268    0.7934    1.3207
   -1.2972   -3.2643   -1.9671    0.7630   -0.1526
    2.2088    5.4597    3.2508    1.3825    0.4301

【讨论】:

  • 非常感谢您的帮助。它有效,但XU 的答案都是[1,...,1][1,...,1],这是不正确的。我认为它来自Unknown = ones(2*r+3, 2*r+1)
  • 我定义了X = Unknown(1, 2*r+1)'; U = Unknown(2*r+2, 4*r+2)';,但我不知道如何定义2*r+1 by 2*r+1矩阵Q作为输入。
  • 你是对的。我问了很多问题,我已经花时间了。您以这种方式定义Unknown = ones(2*r+3, 2*r+1); X = Unknown(1, :)'; U = Unknown(2, :)'; Q = Unknown(3:end, :);,但现在您说删除Unknown = ones(2*r+3, 2*r+1)。我做到了,但它不起作用。我该怎么办?这是我的最后一个问题。非常感谢。
  • 非常感谢您。很抱歉占用您很多时间。非常感谢您的帮助。
  • 您的回答很棒。是的,我当然会。非常感谢。
猜你喜欢
  • 1970-01-01
  • 1970-01-01
  • 2021-12-06
  • 2019-10-02
  • 1970-01-01
  • 1970-01-01
  • 2015-08-21
  • 1970-01-01
  • 2021-12-13
相关资源
最近更新 更多