Can't run an equation in a function that works in a script

Hi,

I've used the following code for an equation in MATLAB and it works perfectly well in a script (cp1 and u1 have already been defined):

syms u4
eqn3 = cp1 == (u1+u4)*(u1.^2-u4.^2)/(2*u1.^3);
U4= double(solve(eqn3,u4));
u4 = real(U4( U4 > 0 ));

However, when I've tried to implement it into a function it doesn't work.

Can anyone help me with this please?

[madhan ravi]

show us your function

[Daniel Harper]

function [y] = test(x)
c1 = 0.5176;
c2 = 116;
c3 = 0.4;
c4 = 5;
c5 = 21;
c6 = 0.0068;
% equation 1 lamdai1
z1 = 1/(x(2) + 0.08*x(1)) - 0.035/((x(1).^3) +1);
y(1)=1/z1;
%equation 2 cp1
y(2) = c1*((c2/y(1))-c3*x(1)-c4)*exp(-c5/y(1))+c6*x(2);
%equation 3  u4 (z3)
syms y(3)
eqn3 = y(2) == (u1+y(3))*(u1.^2-y(3).^2)/(2*u1.^3);
Z3 = double(solve(eqn3,y(3)));
z3 = real(Z3( Z3>0 ));
y(3)=z3
%equation 4
y(4) = 0.5*(u1+y(3));
end

I'm using a genetic algorithm for inputs x(1) and x(2). cp1 and u4 etc. have been equated to y(1) etc.

Answer 1

[Star Strider]

Try this:

function [y] = test(x)
c1 = 0.5176;
c2 = 116;
c3 = 0.4;
c4 = 5;
c5 = 21;
c6 = 0.0068;
% equation 1 lamdai1
z1 = 1/(x(2) + 0.08*x(1)) - 0.035/((x(1).^3) +1);
y = sym('y', [1 4])
y(1)=1/z1;
%equation 2 cp1
y(2) = c1*((c2/y(1))-c3*x(1)-c4)*exp(-c5/y(1))+c6*x(2);
%equation 3  u4 (z3)

u1 = 3;                                                     % <— PROVIDE CORRECT VALUE
eqn3 = y(2) - ((u1+y(3))*(u1.^2-y(3).^2)/(2*u1.^3));
eqn3 = vpa(eqn3);
U4p = sym2poly(eqn3);
Z3 = roots(U4p);
z3 = real(Z3( Z3>0 ));
y(3)=z3;
%equation 4
y(4) = vpa(0.5*(u1+y(3)));
y = double(y);
end

It runs without error for me. You must determine if it produces the correct result.

Note that this cannot be your fitness function, since a fitness function is required to produce a scalar output.

[Daniel Harper]

Thanks it seems to work! However it really slows down the time it takes to run the ga. Is there a way to speed it up?

[Star Strider]

My pleasure!

Speeding it up would require eliminating the Symbolic Math Toolbox calls. I have no idea what you’re doing, so my ability to help you optimize it is significantly restricted.

You can create a polynomial from ‘eqn3’ easily enough by involing the symbolic functions once, then using the result in your code, for example:

syms u1 
y = sym('y', [1 4])
eqn3 = y(2) - ((u1+y(3))*(u1.^2-y(3).^2)/(2*u1.^3))
[eqn3p,trm] = coeffs(eqn3,y(3))

produces:

eqn3p =
[ 1/(2*u1^3), 1/(2*u1^2), -1/(2*u1), y2 - 1/2]

trm =
[ y3^3, y3^2, y3, 1]

(the ‘trm’ output simply shows that the polynomial is sorted correctly).

So ‘Z3’ then becomes:

Z3 = roots([ 1/(2*u1^3), 1/(2*u1^2), -1/(2*u1), y(2) - 1/2]);

So simply use that version of ‘eqn3’ to create the polynomial vector form of it, calculate the roots, and be done with it. It doesn’t appear that any symbolic calls are necessary at all then.

This seems to work:

x = rand(2,1);

c1 = 0.5176;
c2 = 116;
c3 = 0.4;
c4 = 5;
c5 = 21;
c6 = 0.0068;
% equation 1 lamdai1
z1 = 1/(x(2) + 0.08*x(1)) - 0.035/((x(1).^3) +1);
y(1)=1/z1;
%equation 2 cp1
y(2) = c1*((c2/y(1))-c3*x(1)-c4)*exp(-c5/y(1))+c6*x(2);
%equation 3  u4 (z3)
u1 = 3;                                                     % <— PROVIDE CORRECT VALUE
Z3 = roots([ 1/(2*u1^3), 1/(2*u1^2), -1/(2*u1), y(2) - 1/2]);
z3 = real(Z3( Z3>0 ));
y(3)=z3;
%equation 4
y(4) = 0.5*(u1+y(3));

That completely eliminates the symbolic calls. It should be faster.

[Daniel Harper]

Oky great thanks it seems to work much faster! The rest of the script is made up just with other basic equations like eqn. 1.

One more thing I, want to find the maximum optimal values whereas this seems to calculate the mimimum optimal values. How can I correct this?

[Star Strider]

My pleasure.

‘... maximum optimal values ...’

Of what? There is no optimization code in this.

My changes do the same as your original code, just without the symbolic calls.

If you’re talking about your overall ga results, note that all optimization routines (at least that I’m aware of) find the minimum of whatever it is they’re optimizing. If you want to maximize your fitness function, negate its output.

[Daniel Harper]

Yeah sorry I should have been more clearer. This is my full fuction. If I wanted to maximise y(12) how would I do this?

x = rand(4,3,2,1);
rho = 1.225; 
R = 89.15; 
k = 0.04;
xd = 178.3; 
u1 = 11.4; 
c1 = 0.5176;
c2 = 116;
c3 = 0.4;
c4 = 5;
c5 = 21;
c6 = 0.0068;

% equation 1 lamdai1
z1 = 1/(x(2) + 0.08*x(1)) - 0.035/((x(1).^3) +1);
y(1)=1/z1;

%equation 2 cp1
y(2) = c1*((c2/y(1))-c3*x(1)-c4)*exp(-c5/y(1))+c6*x(2);

%equation 3  u4 (z3)
Z3 = roots([ 1/(2*u1^3), 1/(2*u1^2), -1/(2*u1), y(2) - 1/2]);
z3 = real(Z3( Z3>0 ));
y(3)=z3;

%equation 3ii  u2
y(4) = 0.5*(u1+y(3));

%equation 3iii  T
y(5) = rho*pi*(R.^2)*y(4)*(u1-y(3));

%equation 3iv ct
y(6) = (2*y(5))/(rho*pi*(R.^2)*(u1.^2));

%equation 4 vw
y(7) = u1*(1-((R/(k*xd+R)).^2)*(1-sqrt(1-y(6))));

% equation 5i - maximum 10 MW allowed pout1
a = 0.5*y(2)*rho*pi*R.^2*u1.^3;
a(a>10000000)=10000000;
y(8) = a;
% equation 1 lamdai2
z3 = 1/(x(4) + 0.08*x(3)) - 0.035/((x(3).^3) +1);
y(9) = 1/z3;
%equation 2  cp2
y(10) = c1*((c2/y(9))-c3*x(3)-c4)*exp(-c5/y(9))+c6*x(4);

% equation 5ii pout2 
y(11) = 0.5*y(10)*rho*pi*R.^2*y(7).^3;

% equation 6 pout total  
y(12) = y(8) + y(11);

[Star Strider]

Maximize it with respect to what?

The ‘y’ variable is a (1 x 12) vector that is not itself a function of any other variables. What do you intend by ‘maximizing’ it? You can get the largest value (and its index) with the max (link) function.

[Star Strider]

Daniel Harper’s ‘Answer’ moved here:

Maximize it in terms of the inputs x. Like how do I maximize the fitness fuction rather than minimize it? I've seen you can do y = -x for the fuction however that doesn't seem to work

[Star Strider]

If ‘test.m’ is your entire fitness function, it is not going to work. Fitness functions require scalar outputs (at least in my experience), and ‘y’ is a (1 x 12) vector.

Negating the fitness function output is the generally accepted way to maximimize a fitness (or other objective) function. That obviously works if there is a distinct (and preferably bounded) maximum. If your fitness function is for example a function that increases without bound in any direction, there is not going to be a distinct maximum.

So first, plot the function you want to optimize, if possible. (It initially appeared to be a function of two parameters, so that could work.) That will give you some idea of what it’s doing, and if finding a distinct maximum is even possible.

Meanwhile, I have gotten yuouir function to run, and run relatively efficiently. That was the initial objective.