使用 fmincon() 求解一个特殊的非线性程序
Solving a special nonlinear program using fmincon()
我的问题是以下优化问题：
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
其中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 矩阵和 I 是一个 2*r+1 by 1 已知列矩阵。此外，Q 是一个矩阵，它满足
I*I'*(D2')*U=Q*I.
Known matrices: X0, E, D1, D2, Z, P, I
提前致谢。
我的尝试：
Ojective function file name: 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 function file name: 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
Optimization file name: 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=[];
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);
U 和 X 是列向量描述如下 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, :)
In constraint.m there is dimension issue while concatenating
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.m 和 objective_function.m 但请先阅读以了解原因。
Ojective function file name: 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 function file name: 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, :)
Solution
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
我的问题是以下优化问题：
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
其中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 矩阵和 I 是一个 2*r+1 by 1 已知列矩阵。此外，Q 是一个矩阵，它满足
I*I'*(D2')*U=Q*I.
Known matrices:
X0, E, D1, D2, Z, P, I
提前致谢。
我的尝试：
Ojective function file name: 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 function file name: 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
Optimization file name: 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=[];
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);
U 和 X 是列向量描述如下 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, :)
In
constraint.mthere is dimension issue while concatenating
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.m 和 objective_function.m 但请先阅读以了解原因。
Ojective function file name: 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 function file name: 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, :)
Solution
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