【最优化方法】实验四 约束最优化方法的MATLAB实现

实验的目的和要求：通过本次实验使学生较为熟练使用MATLAB软件，并能利用该软件进行约束最优化方法的计算。

实验内容：

１、罚函数法的MATLAB实现

２、可行方向法的MATLAB实现

学习建议：

本次实验就是要通过对一些具体问题的分析进一步熟悉软件的操作并加深对理论知识的理解。

重点和难点：

可行点和辅助函数选取。

一 罚函数法

用罚函数法解min f(x)=(x1-1)^2+x22

​ S.T. g(x)=x2-1>=0

编写下面m文件

fahanshu.m

syms x1 x2 

f=(x1-1)^2+x2^2;

g=x2-1;

x=[x1;x2];

x0=[0;0];

e=0.0001;M=1;

while abs(subs(g,x,x0))>e

  if subs(g,x,x0)<0

​    Q=f+M*g^2;

  else 

​    Q=f;

  end

  x0=xzNewton(Q,x,x0,0.0001);

  M=10*M;

end

x0
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33

运行得：fahanshu

x0 =

1.0000

0.9999

与理论值x=[1;1]很接近。

二 投影梯度法

1.梯度投影法基本原理和步骤

思想：当迭代点是可行域的内点时，将目标函数负梯度作为搜索方向，当迭代点在可行域边界上时，将目标函数负梯度在可行域边界上的投影作为搜索方向。无论何种情况，所构造的方向都是可行下降方向。然后在可行域内沿该方向进行最优一维搜索得到新的迭代点。

MATLAB实现:

2.代码及数值算例:

(1) 程序源代码:

function [ X,FMIN,K ] = tidutouying( f,A,b,x1,x,e )

%  [ X,FMIN,K ] = tidutouying( f,A,b,x1,e ) 梯度投影法

%  f  目标函数

%  A  约束矩阵  b 右端项

%  x1 初始点 x 自由变量

%  e  精度要求

%  X  极小点

%  FMIN 极小值

%  K  迭代次数

%  张超编写与2014/5/3

count=1;

n=length(x1);

tf=jacobian(f,x)';

while 1

  [A1,A2,b1,b2,k]=fenjie(A,b,x1);

  while 1

  M=A1;

  if isempty(M)

​    P=eye(n);

  else

​    P=eye(n)-M'*(M*M')^(-1)*M;

  end

  Pk=-P*subs(tf,x,x1);

  if norm(Pk)<=e

​    if isempty(M)

​      x1;break;

​    else

​      W=(M*M')^(-1)*M*subs(tf,x,x1);

​      u=W;

​      if min(u)<0

​        for i=1:length(u)

​          if u(i)==min(u)

​            j=i;

​          end

​        end

​        A1(j,:)=[];

​      else 

​        x1;break;

​      end

​    end

  else

​    b_=b2-A2*x1;

​    P_=A2*Pk;

​    for i=1:length(P_)

​      if P_(i)<0

​        r(i)=b_(i)/P_(i);

​      else 

​        r(i)=10000;

​      end

​    end

​    rmax=min(r);

​    syms t

​    y=x1+t*Pk;

​    ft(t)=subs(f,x,y);    

​    [r1]=find0618(ft,0,double(rmax),0.00001);

​    x1=x1+r1*Pk;   

​    break;

  end

  end

  count=count+1;

  if isempty(M)    

   break;

  end

  if min(u)>=0

​    break;

  end

end

 X=x1;

 FMIN=subs(f,x,X);

 K=count;

end

 

function [A1,A2,b1,b2,k ] = fenjie( A,b,x )

% 分解起作用约束

A=A;

 b=b;

 x0=x;

 k=0;q=0;

s=size(A);

A1=zeros(s(1),s(2));

 A2=zeros(s(1),s(2));

 b1=zeros(s(1),1);

 b2=zeros(s(1),1);

 for i=1:s(1)

  gi=A(i,:)*x0-b(i);

  if abs(gi)<0.0000001 

​    k=k+1;

​    A1(k,:)=A(i,:);

​    b1(k,1)=b(i);

  else 

​    q=q+1;

​    A2(q,:)=A(i,:);

​    b2(q,1)=b(i);

  end

 end

if k>0

  A1=A1(1:k,:);

  b1=b1(1:k,:);

else

  A1=[];

end

if q>0

  A2=A2(1:q,:);

  b2=b2(1:q,:);

end 

end
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211

(2) 数值算例:

Min f(x)= 2*x1^2+2*x2^2 – 2*x1*x2 – 4*x1 – 6*x2

S.T. –x1 – x2>= –2

–x1 – 5*x2>= –5

x1>=0

x2>=0

初值x0=[0;0]

在matlab command window里输入

syms x1 x2 

f=2*x1^2+2*x2^2-2*x1*x2-4*x1-6*x2;

A=[-1 -1;-1 -5;1 0;0 1];

b=[-2;-5;0;0];

x1=[0;0];

x=[x1;x2];

e=0.01;

[X,FMIN,K]=tidutouying(f,A,b,x1,x,e)

 

X =

  1.1292

  0.7742

FMIN =

  -7.1613

N =

  16
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45