I want to minimize a function, but it gives an error message "Conversion to logical from sym is not possible.". I would like to understand the reason behind it.

At 1st I define te function and then I find the minimum value. The code is as follows:

function q = Bauer_var(v)
syms a m n b r s phi E D ri ro u;
ri = 50;
ro = 200;
E = 210000;
s = 1;
u = 0.3;
D = (E*s^3)/(12*(1-u^2));
a = v(1);
m = v(2);
n = v(3);
b = v(4);
assume(r >= ri);
assume(r <= ro);
z = exp(a*r)*((r-ri)^2)*((r-ro)^2)*cos(n*(tan(b)*log(ri/r) + phi))/(r^m);
Mkrit = (int(int((D*((diff(z,r,2) + (1/r)*diff(z/r) + (1/r^2)*diff(z,phi,2))^2 + 2*(1-u)*(((1/r)*diff(diff(z,r),phi) - (1/r^2)*diff(z,phi))^2 - diff(z,r,2)*((1/r)*diff(z,r) + (1/r^2)*diff(z,phi,2)))))*2,r,[ri ro]),phi,[0 2*pi]))/((-1)*int(int(((1/r^3)*diff(z,r)*diff(z,phi))*r,r,[ri ro]),phi,[0 2*pi]));
q = Mkrit;
end

After defining the function as above, I minimize it as follows:

v = [10, 10, 5, 10];
var_min = fminsearch(@Bauer_var,v);
a = var_min(1)
m = var_min(2)
n = var_min(3)
b = var_min(4)

The above minimisation should be giving me the minimum values, but it always gives the following error:

{{

Conversion to logical from sym is not possible.

Error in fminsearch (line 327)

if fxr < fv(:,1)

Error in BauerVarExec (line 2)

var_min = fminsearch(@Bauer_var,v);

}}

I would really appreciate any help. And I am new to this, so would like to get an explanation which a beginner can unerstand.

Answer 1

[Steven Lord]

The fminsearch function requires the function handle that you pass into it as the first input to return a numeric array, not a symbolic array. Convert the symbolic value Mkrit into double immediately before you return it from your function.

Basically what you're doing is returning sqrt2Symbolic. You need to return sqrt2Double instead.

two = sym(2);
sqrt2Symbolic = sqrt(two)

sqrt2Symbolic = 

sqrt2Double = double(sqrt2Symbolic)

sqrt2Double = 1.4142

[Sarbhanu Saha]

Thank you very much for answering. I am converting Mkrit into double by defininf q = double(Mkrit). Is it the right way? The edited function is written below for reference. Or are you talking of some other way of converting it? Can you please help me understand well how to implement that conversion to double.

function q = Bauer_var(v)
syms a m n b r s phi E D ri ro u;
ri = 50;
ro = 200;
E = 210000;
s = 1;
u = 0.3;
D = (E*s^3)/(12*(1-u^2));
a = v(1);
m = v(2);
n = v(3);
b = v(4);
assume(r >= ri);
assume(r <= ro);
z = exp(a*r)*((r-ri)^2)*((r-ro)^2)*cos(n*(tan(b)*log(ri/r) + phi))/(r^m);
Mkrit = (int(int((D*((diff(z,r,2) + (1/r)*diff(z/r) + (1/r^2)*diff(z,phi,2))^2 + 2*(1-u)*(((1/r)*diff(diff(z,r),phi) - (1/r^2)*diff(z,phi))^2 - diff(z,r,2)*((1/r)*diff(z,r) + (1/r^2)*diff(z,phi,2)))))*2,r,[ri ro]),phi,[0 2*pi]))/((-1)*int(int(((1/r^3)*diff(z,r)*diff(z,phi))*r,r,[ri ro]),phi,[0 2*pi]));
q = double(Mkrit);
end

[Sarbhanu Saha]

If i used the above program, It is taking around 1 day and then is crashing. Can you please explain on how to convert it to double?

[Walter Roberson]

Mkrit has a place with diff(z/r) that should probably be diff(z,r)

[Sarbhanu Saha]

Thank you very much for pointing out an error! I will run again to see if i get any results.

[Walter Roberson]

Reverse the order of integration. The inner integral of the numerator can be symbolically integrated with respect to phi no matter what the v values are.

However, you should probably switch the (new) outer integral to vpaintegral().

It might be worth switching all of the int() to vpaintegral()

[Walter Roberson]

With b = v(4) and the tan(b) in the code, the expressions end up with a fair number of sin(v(4)) and cos(v(4)) terms. You initial v(4) is 10. Could you confirm that you are wanting to explore around 10 radians ?

... Ummm, running the code with your initial values generates an output of about 6E+1703 which is going to overflow as soon as you try to convert it to double precision...

[Walter Roberson]

Tests at the moment are showing hints of V1=0, V2=0, V3=0.35621657788, V4=0, but I have to do more checking.

[Sarbhanu Saha]

I tried to check again the equation and corrected all the errors. The updated code is as follows:

function q = Bauer_var(v)
syms a m n b r s phi E D ri ro u;
ri = 50;
ro = 200;
E = 210000;
s = 1;
u = 0.3;
D = (E*s^3)/(12*(1-u^2));
a = v(1);
m = v(2);
n = v(3);
b = v(4);
assume(r >= ri);
assume(r <= ro);
z = exp(a*r)*((r-ri)^2)*((r-ro)^2)*cos(n*(tan(b)*log(ri/r) + phi))/(r^m);
Mkrit = (pi*vpaintegral(vpaintegral((D*((diff(z,r,2) + (1/r)*diff(z,r) + (1/r^2)*diff(z,phi,2))^2 + 2*(1-u)*(((1/r)*diff(diff(z,r),phi) - (1/r^2)*diff(z,phi))^2 - diff(z,r,2)*((1/r)*diff(z,r) + (1/r^2)*diff(z,phi,2)))))*r,r,[ri ro]),phi,[0 2*pi]))/((-1)*vpaintegral(vpaintegral(((1/r^3)*diff(z,r)*diff(z,phi))*r,r,[ri ro]),phi,[0 2*pi]));
q = Mkrit;
end

And then minimization is as follows:

v = [10, 10, 5, 10];
var_min = fminsearch(@Bauer_var,v);
a = var_min(1)
m = var_min(2)
n = var_min(3)
b = var_min(4)

[Sarbhanu Saha]

@Walter Roberson

Thank you very much for helping. Actually the value 10 is really a guess value. And we can put any value to be true.

The overfowing is maybe becuse of the difference in power of r in numerator and denominator, which i corrected after editing the equation of Mkrit

Can you please explain how you got the value of v1, v2, v3 and v4... it always ends with an error for me?

And with the new codes using vpaintegral, is it necessary to convert the Mkrit from symbolic fn to Double?

[Walter Roberson]

Here is my current version of the code.

If you study the structure carefully you will notice that for the denominator, I reversed the order of integration. The inner part of the denominator has a closed form integral with respect to phi.

In your previous version of the expressions, the numerator also had a closed form with respect to phi.

This code uses a relatively new option to postpone integration.

I just noticed you are using R2017b; that option probably does not exist as yet. For your version, instead of being able to use int() with hold, and decide afterwards how to turn the formula into an executable function, it might be necessary to convert the last three int() directly into vpaintegral() -- if you leave them as int() without the 'hold' option, then MATLAB will think about them for about 3/4 of an hour before deciding that it does not know how to integrate them :(

syms V [1 4]
tic
formula = Bauer_var(V);
toc

%fun = matlabFunction(formula, 'vars', {V});
fcn(V) = vpa(formula);
fun = @(V) double(fcn(V(1),V(2),V(3),V(4)));

options = optimset(@fminsearch);
options = optimset(options, 'Display', 'iter');

V0 = [0.1, 0.2, 0.3, 0.4];
[bestv, fval] = fminsearch(fun, V0, options);

function q = Bauer_var(v)
syms a m n b r s phi E D ri ro u;
Q = @(v) sym(v);
ri = Q(50);
ro = Q(200);
E = Q(210000);
s = Q(1);
u = Q(3)./10;
D = (E.*s.^3)./(Q(12).*(1-u.^2));
a = v(1);
m = v(2);
n = v(3);
b = v(4);
Pi = sym(pi);
assume(r >= ri);
assume(r <= ro);
z = exp(a.*r).*((r-ri).^2).*((r-ro).^2).*cos(n.*(tan(b).*log(ri./r) + phi))./(r.^m);
dzr = diff(z,r);
dzrr = diff(dzr,r);
dzp = diff(z,phi);
dzpp = diff(dzp,phi);
dzrp = diff(dzr,phi);
q2inside = ((1./r.^3).*dzr.*dzp).*r;
q2inint = int(q2inside,phi,[0 2.*Pi]);          %this one can be integrated
q2 = (-1).*int(q2inint,r,ri,ro, 'Hold', true);

q1inside = (D.*((dzrr + (1./r).*dzr + (1./r.^2).*dzpp).^2 + 2.*(1-u).*(((1./r).*dzrp - (1./r.^2).*dzp).^2 - dzrr.*((1./r).*dzr + (1./r.^2).*dzpp)))).*r;
q1inint = int(q1inside,r,ri, ro, 'Hold', true);
q1 = Pi.*int(q1inint,phi,0, 2.*Pi, 'Hold', true);
Mkrit = q1./q2;
q = Mkrit;
end

[Sarbhanu Saha]

Thank you very much for putting so much effort into it. I wanted to know few things, can you help me clear it.....

What is the use of the Q() function that you applied on constans?

What is the benifit of using .* and ./ rather than directly using * or /?

Does these two effect the speed of the operation anyway?

[Walter Roberson]

Q() is a little helper function to make it easier to convert between symbolic numbers and floating point numbers -- and to go back to floating point easily if desired.

The code for it is given right there at the top

Q = @(v) sym(v);

so Q(50) is sym(50) which produces the symbolic number 50 from the floating point number 50. The name Q was chosen because ℚis the mathematical symbol for the domain of rational numbers, and most of the time sym() converts floating point values into rational values. (It can detect simple multiples of pi, and can detect square roots of rationals up to about 5 digits, so it does not always produce rationals.)

When you are working with symbolic forms and with floating point numbers, it is often best to convert the floating point numbers to symbolic form early. For example, if you had a term

syms t
sin(2/456789*pi*t)

ans = 

then MATLAB would interpret 2 and pi as floating point values and calculate the result, and would divide the result by 456789 in floating point. Then, it would see the symbolic t and decide that it needed to convert the double precision number. Using the default conversion algorithm it would convert it to the rational number 8119586860750777/590295810358705651712 .

If, on the other hand, you have

sin(2/456789*sym(pi)*t)

ans = 

sin(sym(2)/456789*pi*t)

ans = 

It is easier to be mindless and convert all the floating point numbers as they are introduced, instead of analyzing whether a particular section of code will happen to build the symbolic numbers "nicely"... especially since if you are consistent, then if you edit the code later you do not need to think about whether there is anything magic about the order of constants in the code. sym() them all to get predictable values.

Using .* and ./ instead of * and / often does not matter in symbolic code, as much of the time symbolic expressions turn out to be scalar. Until R2021a, symbolic variable names always stood in for scalars. But.. people using symbolics tend to lose track of the fact that they introduced a numeric vector and so suddenly they are working with vectors of expressions and having odd results when they multiply or divide...

What is the benifit of using .* and ./ rather than directly using * or /?

In MATLAB, it is the other way around. Element-by-element multiplication is the .* operator and * is algebraic matrix multiplication (inner product). 3*5 is not a matter of "directly" multiplying 3 and 5: 3*5 has to check the size of variables and notice that at least one of them is scalar, and call over to the element-by-element multiplication in that case: the .* operator is the "real" operator in that case, with the * operator being effectively a permitted abbreviation.

So... when I am editing code for people, I routinely change to .* and ./ and .^ (except sometimes for simple scalar constants.) It tends to fix problems people have overlooked, and makes it easier to extend the code to vectors in the future.

In the case of your code, I had the additional motivation that I was making it easier to convert the code to use numeric vectors, with the intent of seeing if I could speed it up by vectorization. I did not take the final step of converting the int() and nested int() into integral() and integral2() calls, but the infrastructure is better prepared for that.

[Walter Roberson]

I am running a minimization based upon the code, using a starting point of [0.1, 0.2, 0.3, 0.4] . After 90 iterations (155 function calls) it descended from about 4E+7 "cost" down to about -9e9, and is currently around -3e14 . However, at the moment I do not know what range it is exploring.

[Walter Roberson]

It is now taking about 3 hours per iteration (with it doing two function calls per iteration.) I guess it must be having trouble doing the integrations, suggesting that either it is approaching discontinuities or else that it has entered an area that is changing rapidly.

[Sarbhanu Saha]

@Walter Roberson You explained it really in a simple way and so i could understand both Q() and .* & ./. Thank you very much for giving time in explaining me.

I have run the code that was before your correction and modification, and it has been running since 4th July evening, fortunately without errors. After this, I will run your modified code using vpaintegral replacing the Hold function, since it is R2017b.

The word "Cost" that you are using, is it the value of some kind of error or deviation from the real value or deviation between some calculated values? And when the algorithm is programmed to stop, like at what value of "Cost"? Can you please explain? And can we decide that final value of "Cost" to stop the algorithm?

[Walter Roberson]

In minimization, "cost" refers to the value being minimized. In particular the function being minimized is referred to as the "cost function" . It is easier to say things like "The cost function will initially be evaluated (N+1) times for N variables" than it is to keep having to say "The function that you passed as the first parameter will initially be evaluated (N+1) times for N parameters". "cost function" describes the function by role, and so is much better suited for discussion in theory papers instead of referring to the implementation detail that the function happens to be the first parameter to a particular MATLAB library routine.

Evaluation stops under several different circumstances:

• configured limit on the maximum number of function evaluations is reached
• configured limit on the number of iterations is reached
• for some algorithms, the first time the function returns nan or infinity or complex numbers. Other algorithms can continue trying
• the algorithm reaches a point where all of its strategies for moving result in values that are the same (fallen into a flat area)
• the algorithm reaches a point where all of its strategies for moving result in values that are the same or higher (global minimum or steep enough local minimum)

[Walter Roberson]

After about 2 days, about 215 function calls, it gets to the point where vpaintegral() fails. Unfortunately with the current code, that causes minimization to end, losing all progress.

To get around that, the code

fun = @(V) double(fcn(V(1),V(2),V(3),V(4)));

needs to be replaced with a non-anonymous function, something along the lines of

function cost = failable_fcn(V, fcn)
   cost = 1e50;
   cost_sym = [];
   try
       cost_sym = fcn(V);
       cost = double(cost_sym);   %might fail for nonconvergence
   catch ME
       if isempty(cost_sym)
          fprintf('fcn failed completely at ');
          disp(V);
       else
          fprintf('fcn worked, double() failed, at ');
          disp(V);
          disp(cost_sym);
       end
   end
end

which would be used like

fcn(V) = vpa(formula);
fun = @(V) failable_fcn(V, fcn);

[Sarbhanu Saha]

@Walter Roberson

Thank you very much for explaining the "Cost" value.

I have been trying the use of failable_fcn, but it always gives the same value what i enter as the initial value whatsoever.

There is one condition which i can apply on "n", which is in the coding is v(3). It was supposed to come as an integer between 1 and 10. Can that help in any way?

[Walter Roberson]

I created a version for testing, in which I converted to using vpaintegral() at each stage. I invoke it inside of ga() since ga() is able to handle integer constraints while it minimizes.

Unfortunately, the version with vpaintegral() is wretchedly slow. It is able to create the symbolic function in a quite decent time, in terms of a series of vpaintegral() calls, but the resulting function is too slow to get anywhere. I have had it executing for about 3 hours but it has yet to finish a single evaluation.

Yes, it is true that in the previous version, it eventually slowed down to about that speed, but this is the very first evaluation, which only took about 3 minutes in the previous version.

I would firmly recommend that you upgrade to a modern MATLAB; the methods needed for your older release are not feasible.

[Sarbhanu Saha]

@Walter Roberson

I have been working on it. And i now found the minimisation values of the function by getting the absolute value of Mkrit by abs(Mkrit). From the help of my PhD guides, I also changed the initial vues for the iteration. Here is the execution:

syms V [1 4]
formula = Bauer_var(V);
fcn(V) = vpa(formula);
fun = @(V) double(fcn(V(1),V(2),V(3)));
options = optimset(@fminsearch);
options = optimset(options, 'Display', 'iter');
V0 = [0.005, 2, 1, pi/4];
[bestv, fval] = fminsearch(fun, V0, options);
a = bestv(1)
m = bestv(2)
n = bestv(3)
b = bestv(4)

Below is the function, which will be minimised:

function q = Bauer_var(v)
tic
syms a m n b r s phi E D ri ro u;
Q = @(v) sym(v);
ri = Q(10);
ro = Q(80);
E = Q(210000);
s = Q(1);
u = Q(3)./10;
D = (E.*s.^3)./(Q(12).*(1-u.^2));
a = v(1);
m = v(2);
n = v(3);
b = v(4);
Pi = sym(pi);
assume(r >= ri);
assume(r <= ro);
z = exp(a.*r).*((r-ri).^2).*((r-ro).^2).*cos(n.*(tan(b).*log(ri./r) + phi))./(r.^m);
dzr = diff(z,r);
dzrr = diff(dzr,r);
dzp = diff(z,phi);
dzpp = diff(dzp,phi);
dzrp = diff(dzr,phi);
q2inside = ((1./r.^3).*dzr.*dzp).*r;
q2inint = vpaintegral(q2inside,phi,[0 2.*Pi]);          
q2 = (-1).*vpaintegral(q2inint,r,[ri ro]);
q1inside = (D.*((dzrr + (1./r).*dzr + (1./r.^2).*dzpp).^2 + 2.*(1-u).*(((1./r).*dzrp - (1./r.^2).*dzp).^2 - dzrr.*((1./r).*dzr + (1./r.^2).*dzpp)))).*r;
q1inint = vpaintegral(q1inside,r,[ri ro]);
q1 = Pi.*vpaintegral(q1inint,phi,[0, 2.*Pi]);
Mkrit = q1./q2;
q = abs(Mkrit);
toc
end