1 | function diagnostic = callbmialt(F,h,options) |
---|
2 | % |
---|
3 | % EXTREMELY naive implementation of a local BMI solver |
---|
4 | % using alternating directions. |
---|
5 | % |
---|
6 | % The default behaviour of the solver can be |
---|
7 | % altered be using the field bmialt in sdpsettings. |
---|
8 | % |
---|
9 | |
---|
10 | % Recover all used variables |
---|
11 | |
---|
12 | x_vars = depends(F); |
---|
13 | all_vars = getvariables(F); |
---|
14 | |
---|
15 | % Nonlinear variables |
---|
16 | nonlin_variables = setdiff(all_vars,x_vars); |
---|
17 | |
---|
18 | x_nonlin = recover(nonlin_variables); |
---|
19 | |
---|
20 | set1_vars = []; |
---|
21 | set2_vars = []; |
---|
22 | for i = 1:length(x_nonlin); |
---|
23 | xi_vars = depends(x_nonlin(i)); |
---|
24 | if length(xi_vars)~=2 |
---|
25 | error('Not bilinear') |
---|
26 | end |
---|
27 | if isempty(set1_vars) |
---|
28 | set1_vars = xi_vars(1); |
---|
29 | set2_vars = xi_vars(2); |
---|
30 | else |
---|
31 | % index = [ismember(xi_vars(1),set1_vars) ismember(xi_vars(2),set1_vars) ismember(xi_vars(1),set2_vars) ismember(xi_vars(2),set2_vars)]; |
---|
32 | |
---|
33 | if ismember(xi_vars(1),set1_vars) |
---|
34 | set2_vars = [set2_vars xi_vars(2)]; |
---|
35 | elseif ismember(xi_vars(1),set2_vars) |
---|
36 | set1_vars = [set1_vars xi_vars(2)]; |
---|
37 | elseif ismember(xi_vars(2),set1_vars) |
---|
38 | set2_vars = [set2_vars xi_vars(1)]; |
---|
39 | elseif ismember(xi_vars(2),set2_vars) |
---|
40 | set1_vars = [set1_vars xi_vars(1)]; |
---|
41 | end |
---|
42 | end |
---|
43 | end |
---|
44 | |
---|
45 | |
---|
46 | set1_vars = unique(set1_vars); |
---|
47 | set2_vars = unique(set2_vars); |
---|
48 | |
---|
49 | x = recover(x_vars); |
---|
50 | x1 = recover(set1_vars); |
---|
51 | x2 = recover(set2_vars); |
---|
52 | |
---|
53 | % Set up for local solver |
---|
54 | verbose = options.verbose; |
---|
55 | options.verbose = max(options.verbose-1,0); |
---|
56 | options.solver = options.bmilin.solver; |
---|
57 | |
---|
58 | % Initial values hopefully supplied |
---|
59 | if options.usex0 |
---|
60 | % Initialize to 0 if not initialized |
---|
61 | not_initial = isnan(double(x)); |
---|
62 | if any(not_initial) |
---|
63 | setsdpvar(x(find(not_initial)),repmat(0,length(find(not_initial)),1)); |
---|
64 | end |
---|
65 | else |
---|
66 | % No initials, set to zero |
---|
67 | setsdpvar(x,repmat(0,length(x),1)); |
---|
68 | F_linear = F(find(is(F,'linear'))); |
---|
69 | % Find some non-zero by solving for the linear constraints |
---|
70 | if length(F_linear) > 0 |
---|
71 | solvesdp(F_linear,linearize(h)+x'*x,options); |
---|
72 | end |
---|
73 | end |
---|
74 | |
---|
75 | |
---|
76 | % Outer linearization loop |
---|
77 | goon = 1; |
---|
78 | outer_iter = 0; |
---|
79 | while goon |
---|
80 | |
---|
81 | % Save old iterates and old objective function |
---|
82 | x0 = double(x); |
---|
83 | h0 = double(h); |
---|
84 | |
---|
85 | Flin0 = replace(F,x1,double(x1)); |
---|
86 | obj1 = replace(h,x1,double(x1)); |
---|
87 | |
---|
88 | Flin = set([]); |
---|
89 | for i = 1:length(Flin0) |
---|
90 | if isa(sdpvar(Flin0(i)),'sdpvar') |
---|
91 | Flin = Flin + Flin0(i); |
---|
92 | end |
---|
93 | end |
---|
94 | |
---|
95 | % Solve linearized problem |
---|
96 | solvesdp(Flin,obj1,options); |
---|
97 | |
---|
98 | if verbose > 0 |
---|
99 | disp(sprintf('#%d cost : %6.3f ',outer_iter,double(h))); |
---|
100 | end |
---|
101 | |
---|
102 | Flin0 = replace(F,x2,double(x2)); |
---|
103 | obj1 = replace(h,x2,double(x2)); |
---|
104 | |
---|
105 | Flin = set([]); |
---|
106 | for i = 1:length(Flin0) |
---|
107 | if isa(sdpvar(Flin0(i)),'sdpvar') |
---|
108 | Flin = Flin + Flin0(i); |
---|
109 | end |
---|
110 | end |
---|
111 | |
---|
112 | |
---|
113 | % Solve linearized problem |
---|
114 | solvesdp(Flin,obj1,options); |
---|
115 | |
---|
116 | outer_iter = outer_iter + 1; |
---|
117 | if verbose > 0 |
---|
118 | disp(sprintf('#%d cost : %6.3f ',outer_iter,double(h))); |
---|
119 | end |
---|
120 | goon = (outer_iter < options.bmilin.maxiter); |
---|
121 | end |
---|
122 | |
---|
123 | diagnostic.solvertime = 0; |
---|
124 | diagnostic.info = yalmiperror(0,'BMILIN'); |
---|
125 | diagnostic.problem = 0; |
---|