Solution does not change with different inputs?

1

I am replicating the following equation 7 from Daria Nesterovich Anderson et al 2018 J. Neural Eng. 15 026005:

Screen Shot 2022-03-01 at 1.45.49 PM
Screen Shot 2022-03-01 at 1.45.49 PM776×476 70.6 KB

When I implement my code, my optimization is always the same no matter what I use as my target region, omega. Any thoughts on where I’ve gone wrong? Thank you.

My code is below:

size_omega = sum(sum(sum(omega)));
Integral_omega = zeros(1,8);


for i=1:8
    if (i==1)
        
        [contact1] = voltagebox(contact1);
        [all_derivatives] = Part2ndDeriv(contact1);
        dy2 = all_derivatives{2,1};
        dx2 = all_derivatives{2,2};
        dxy = all_derivatives{2,4};
        dyx = all_derivatives{2,7};
        dz2 = all_derivatives{2,3};
        dxz = all_derivatives{2,5};
        dzx = all_derivatives{2,8};
        dyz = all_derivatives{2,6};
        dzy = all_derivatives{2,9};
    elseif (i==2)
        
        [contact2] = voltagebox(contact2);
        [all_derivatives] = Part2ndDeriv(contact2);
        dy2 = all_derivatives{2,1};
        dx2 = all_derivatives{2,2};
        dxy = all_derivatives{2,4};
        dyx = all_derivatives{2,7};
        dz2 = all_derivatives{2,3};
        dxz = all_derivatives{2,5};
        dzx = all_derivatives{2,8};
        dyz = all_derivatives{2,6};
        dzy = all_derivatives{2,9};
    
    elseif (i==3)
        
        [contact3] = voltagebox(contact3);
        [all_derivatives] = Part2ndDeriv(contact3);
        dy2 = all_derivatives{2,1};
        dx2 = all_derivatives{2,2};
        dxy = all_derivatives{2,4};
        dyx = all_derivatives{2,7};
        dz2 = all_derivatives{2,3};
        dxz = all_derivatives{2,5};
        dzx = all_derivatives{2,8};
        dyz = all_derivatives{2,6};
        dzy = all_derivatives{2,9};
    
    elseif (i==4)
        
        [contact4] = voltagebox(contact4);
        [all_derivatives] = Part2ndDeriv(contact4);
        dy2 = all_derivatives{2,1};
        dx2 = all_derivatives{2,2};
        dxy = all_derivatives{2,4};
        dyx = all_derivatives{2,7};
        dz2 = all_derivatives{2,3};
        dxz = all_derivatives{2,5};
        dzx = all_derivatives{2,8};
        dyz = all_derivatives{2,6};
        dzy = all_derivatives{2,9};
    
    elseif (i==5)
        
        [contact5] = voltagebox(contact5);
        [all_derivatives] = Part2ndDeriv(contact5);
        dy2 = all_derivatives{2,1};
        dx2 = all_derivatives{2,2};
        dxy = all_derivatives{2,4};
        dyx = all_derivatives{2,7};
        dz2 = all_derivatives{2,3};
        dxz = all_derivatives{2,5};
        dzx = all_derivatives{2,8};
        dyz = all_derivatives{2,6};
        dzy = all_derivatives{2,9};

    elseif (i==6)
        
        [contact6] = voltagebox(contact6);
        [all_derivatives] = Part2ndDeriv(contact6);
        dy2 = all_derivatives{2,1};
        dx2 = all_derivatives{2,2};
        dxy = all_derivatives{2,4};
        dyx = all_derivatives{2,7};
        dz2 = all_derivatives{2,3};
        dxz = all_derivatives{2,5};
        dzx = all_derivatives{2,8};
        dyz = all_derivatives{2,6};
        dzy = all_derivatives{2,9};     

    elseif (i==7)
        
        [contact7] = voltagebox(contact7);
        [all_derivatives] = Part2ndDeriv(contact7);
        dy2 = all_derivatives{2,1};
        dx2 = all_derivatives{2,2};
        dxy = all_derivatives{2,4};
        dyx = all_derivatives{2,7};
        dz2 = all_derivatives{2,3};
        dxz = all_derivatives{2,5};
        dzx = all_derivatives{2,8};
        dyz = all_derivatives{2,6};
        dzy = all_derivatives{2,9};

elseif (i==8)
        
        [contact8] = voltagebox(contact8);
        [all_derivatives] = Part2ndDeriv(contact8);
        dy2 = all_derivatives{2,1};
        dx2 = all_derivatives{2,2};
        dxy = all_derivatives{2,4};
        dyx = all_derivatives{2,7};
        dz2 = all_derivatives{2,3};
        dxz = all_derivatives{2,5};
        dzx = all_derivatives{2,8};
        dyz = all_derivatives{2,6};
        dzy = all_derivatives{2,9};        

    end
        for x=1:199
            for y=1:199
                for z=1:199
                    if (omega(x,y,z)==1) %if (x,y,z) is a point in the ROI omega:
                        % Find the hessian at (x,y,z):
                        hess = [dx2(x,y,z),dxy(x,y,z),dxz(x,y,z); dyx(x,y,z),dy2(x,y,z),dyz(x,y,z); dzx(x,y,z),dzy(x,y,z),dz2(x,y,z)];
                            a=(isnan(hess) | isinf(hess));
                            if (sum(sum(a)>0))
                                [m n] = find(a==1);
                                aa = length(m);
                                for b=1:aa
                                    hess(m(b),n(b)) = 0;
                                end
                            end
        
                        % Determine the eigenvalues of the hessian:
                        [~, val] = eig(hess);
        
                        % Sum the eigenvalues:
                        lambda_sum = val(1)+val(2)+val(3);
        
                        Integral_omega(i) = lambda_sum + Integral_omega(i);
        
                    else 
                        Integral_omega(i) = Integral_omega(i) + 0;
                    end
                end
            end
        end
end

cvx_begin
    variables c(i) 
    maximize sum(c(i)/size_omega * Integral_omega)
    subject to
    c<=4
cvx_end
3

I have no idea what you’re doing, or even what you ran, because you don’t show cvx_begin to cvx_end as being in a for loop. And you have not provided reproducible code.

If you run this after the preceding for loop, i will have the value 8. You are declaring a variable c having i elements, but only the last of these elements contributes to the objective function. Perhaps you want maximize(c' * Integral_omega/size_omega) or equivalently maximize(sum(c/size_omega .* Integral_omega))

I haven’t looked at what you’ve done, so there might be more wrong than this. I don’t know what you changed when omega is changed, so I can’t comment on why “the optimization” (whatever you mean by that) is always the same. If Integral_omega consists of only nonnegative elements (and you make the change I suggested), then it will always be optimal for all elements of c to be 4 (per your CVX program - I have no idea what should happen in the problem you want to solve). With the program as is, only c(i) is used in the objective function, not all elements of c, and optimal value of c(i) would always be 4.

And if, in my suggested modified version, any element of Integral_omega is negative, then the optimization problem is unbounded, because the corresponding element of c can be chosen as -\infty. With your program as written, the problem is unbounded if Integral_omega(i) < 0. So your optimization problem formulation appears to be not very good, and needs to be re-thought.