-
Notifications
You must be signed in to change notification settings - Fork 1
/
mmasub.m
159 lines (149 loc) · 5.23 KB
/
mmasub.m
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
%-------------------------------------------------------
% This is the file mmasub.m
%
function [xmma,ymma,zmma,lam,xsi,eta,mu,zet,s,low,upp] = ...
mmasub(m,n,iter,xval,xmin,xmax,xold1,xold2, ...
f0val,df0dx,fval,dfdx,low,upp,a0,a,c,d);
%
% Version September 2007 (and a small change August 2008)
%
% Krister Svanberg <[email protected]>
% Department of Mathematics, SE-10044 Stockholm, Sweden.
%
% This function mmasub performs one MMA-iteration, aimed at
% solving the nonlinear programming problem:
%
% Minimize f_0(x) + a_0*z + sum( c_i*y_i + 0.5*d_i*(y_i)^2 )
% subject to f_i(x) - a_i*z - y_i <= 0, i = 1,...,m
% xmin_j <= x_j <= xmax_j, j = 1,...,n
% z >= 0, y_i >= 0, i = 1,...,m
%*** INPUT:
%
% m = The number of general constraints.
% n = The number of variables x_j.
% iter = Current iteration number ( =1 the first time mmasub is called).
% xval = Column vector with the current values of the variables x_j.
% xmin = Column vector with the lower bounds for the variables x_j.
% xmax = Column vector with the upper bounds for the variables x_j.
% xold1 = xval, one iteration ago (provided that iter>1).
% xold2 = xval, two iterations ago (provided that iter>2).
% f0val = The value of the objective function f_0 at xval.
% df0dx = Column vector with the derivatives of the objective function
% f_0 with respect to the variables x_j, calculated at xval.
% fval = Column vector with the values of the constraint functions f_i,
% calculated at xval.
% dfdx = (m x n)-matrix with the derivatives of the constraint functions
% f_i with respect to the variables x_j, calculated at xval.
% dfdx(i,j) = the derivative of f_i with respect to x_j.% low = Column vector with the lower asymptotes from the previous
% iteration (provided that iter>1).
% upp = Column vector with the upper asymptotes from the previous
% iteration (provided that iter>1).
% a0 = The constants a_0 in the term a_0*z.
% a = Column vector with the constants a_i in the terms a_i*z.
% c = Column vector with the constants c_i in the terms c_i*y_i.
% d = Column vector with the constants d_i in the terms 0.5*d_i*(y_i)^2.
%
%*** OUTPUT:
%
% xmma = Column vector with the optimal values of the variables x_j
% in the current MMA subproblem.
% ymma = Column vector with the optimal values of the variables y_i
% in the current MMA subproblem.
% zmma = Scalar with the optimal value of the variable z
% in the current MMA subproblem.
% lam = Lagrange multipliers for the m general MMA constraints.
% xsi = Lagrange multipliers for the n constraints alfa_j - x_j <= 0.
% eta = Lagrange multipliers for the n constraints x_j - beta_j <= 0.
% mu = Lagrange multipliers for the m constraints -y_i <= 0.
% zet = Lagrange multiplier for the single constraint -z <= 0.
% s = Slack variables for the m general MMA constraints.
% low = Column vector with the lower asymptotes, calculated and used
% in the current MMA subproblem.
% upp = Column vector with the upper asymptotes, calculated and used
% in the current MMA subproblem.
%
%epsimin = sqrt(m+n)*10^(-9);
epsimin = 10^(-7);
raa0 = 0.00001;
move = 1.0;
%move = 0.5;
albefa = 0.1;
%albefa = 0.05;
asyinit = 0.5;
asyinit = 0.125;
asyincr = 1.2;
asyincr = 1.1;
asydecr = 0.7;
asydecr = 0.85;
eeen = ones(n,1);
eeem = ones(m,1);
zeron = zeros(n,1);
% Calculation of the asymptotes low and upp :
if iter < 2.5
low = xval - asyinit*(xmax-xmin);
upp = xval + asyinit*(xmax-xmin);
else
zzz = (xval-xold1).*(xold1-xold2);
factor = eeen;
factor(find(zzz > 0)) = asyincr;
factor(find(zzz < 0)) = asydecr;
low = xval - factor.*(xold1 - low);
upp = xval + factor.*(upp - xold1);
lowmin = xval - 10*(xmax-xmin);
lowmax = xval - 0.01*(xmax-xmin);
uppmin = xval + 0.01*(xmax-xmin);
uppmax = xval + 10*(xmax-xmin);
low = max(low,lowmin);
low = min(low,lowmax);
upp = min(upp,uppmax);
upp = max(upp,uppmin);
end
% Calculation of the bounds alfa and beta :
zzz1 = low + albefa*(xval-low);
zzz2 = xval - move*(xmax-xmin);
zzz = max(zzz1,zzz2);
alfa = max(zzz,xmin);
zzz1 = upp - albefa*(upp-xval);
zzz2 = xval + move*(xmax-xmin);
zzz = min(zzz1,zzz2);
beta = min(zzz,xmax);
% Calculations of p0, q0, P, Q and b.
xmami = xmax-xmin;
xmamieps = 0.00001*eeen;
xmami = max(xmami,xmamieps);
xmamiinv = eeen./xmami;
ux1 = upp-xval;
ux2 = ux1.*ux1;
xl1 = xval-low;
xl2 = xl1.*xl1;
uxinv = eeen./ux1;
xlinv = eeen./xl1;
%
p0 = zeron;
q0 = zeron;
p0 = max(df0dx,0);
q0 = max(-df0dx,0);
%p0(find(df0dx > 0)) = df0dx(find(df0dx > 0));
%q0(find(df0dx < 0)) = -df0dx(find(df0dx < 0));
pq0 = 0.001*(p0 + q0) + raa0*xmamiinv;
p0 = p0 + pq0;
q0 = q0 + pq0;
p0 = p0.*ux2;
q0 = q0.*xl2;
%
P = sparse(m,n);
Q = sparse(m,n);
P = max(dfdx,0);
Q = max(-dfdx,0);
%P(find(dfdx > 0)) = dfdx(find(dfdx > 0));
%Q(find(dfdx < 0)) = -dfdx(find(dfdx < 0));
PQ = 0.001*(P + Q) + raa0*eeem*xmamiinv';
P = P + PQ;
Q = Q + PQ;
P = P * spdiags(ux2,0,n,n);
Q = Q * spdiags(xl2,0,n,n);
b = P*uxinv + Q*xlinv - fval ;
%
%%% Solving the subproblem by a primal-dual Newton method
[xmma,ymma,zmma,lam,xsi,eta,mu,zet,s] = ...
subsolv(m,n,epsimin,low,upp,alfa,beta,p0,q0,P,Q,a0,a,b,c,d);