[37] | 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; |
---|