-
Notifications
You must be signed in to change notification settings - Fork 0
/
generateAnalyticalOPCARealization.m
59 lines (59 loc) · 2.34 KB
/
generateAnalyticalOPCARealization.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
function [x, S] = generateAnalyticalOPCARealization(U, Sig, reducedDim, randVector, xm, gamma)
% Authors: H. X. Vo and L. J. Durlorfky
%% This function is not written for efficiency but rather than for easy to understand
a = U(:,1:reducedDim) * diag(Sig(1:reducedDim)) * randVector + xm; % PCA solution
x = zeros(length(a), 1);
mu = zeros(length(a), 1); % Lagrange multiplier
eta = zeros(length(a), 1); % Lagrange multiplier
LagrangianGrad = zeros(1, length(a));
S = zeros(length(a), length(randVector)); % sensitivity matrix: S = dx/dxi
for i=1:length(a)
if gamma > 1 % this scenario is not used as we will get a concave combined objective function
f0 = (1-gamma) * 0 - 2 * (a(i) - gamma/2) * 0;
f1 = (1-gamma) * 1 - 2 * (a(i) - gamma/2) * 1;
if f1 < f0
x(i) = 1;
mu(i) = 0;
eta(i) = 2 * a(i) + gamma - 2;
elseif f0 <= f1
x(i) = 0;
mu(i) = gamma - 2 * a(i);
eta(i) = 0;
end
elseif gamma == 1 % this scenario is not used as we will get a linear combined objective function
f0 = -2 * (a(i) - gamma/2) * 0;
f1 = -2 * (a(i) - gamma/2) * 1;
if f1 < f0
x(i) = 1;
mu(i) = 0;
eta(i) = 2 * a(i) + gamma - 2;
elseif f0 <= f1
x(i) = 0;
mu(i) = gamma - 2 * a(i);
eta(i) = 0;
end
else % we use this scenario because we want a convex combined objective function
xS = (a(i) - gamma/2) / (1 - gamma);
if (xS < 1) && (xS > 0)
x(i) = xS;
mu(i) = 0;
eta(i) = 0;
elseif (xS >= 1)
x(i) = 1;
mu(i) = 0;
eta(i) = 2 * a(i) + gamma - 2;
else
x(i) = 0;
mu(i) = gamma - 2 * a(i);
eta(i) = 0;
end
for j=1:length(randVector)
da_idxi_j = Sig(j)*U(i, j);
S(i, j) = (2*da_idxi_j * x(i) * (1 - x(i))) /...
(mu(i) + 2 * x(i) + eta(i) * x(i) - 2 * gamma * x(i) - mu(i) * x(i) + 2 * gamma * x(i)^2 - 2 * x(i)^2); % sensitivity matrix dx/dxi
end
end
LagrangianGrad(i) = 2 * (1 - gamma) * x(i) - 2 * (a(i) - gamma/2) - mu(i) + eta(i);
end
fprintf('L2 Norm of Lagrangian Gradient = %f\n', norm(LagrangianGrad,2));
end