Evaulate Piecewise MATLAB function

Hi, I have a function with a set of equations which needs to be solved for unknowns x. However, this function consists of piecewise function of matlab. The code looks something like

function dx = myFunction(x)
  a1 = x(1);
  a2 = x(2);

    c = piecewise(a2>0,a1^2,0)

     dx(1) = a1*a2 + a2*c + a1*c;
     dx(2) = a2*c;

end

When I try to use fsolve on this function, MATLAB gives an error message saying the following: Undefined function 'piecewise' for input arguments of type 'logical'

I was wondering if there is a way I can make MATLAB evaluate the piecewise expression for each iteration of fsolve.

Thanks, Rupa

Answer 1

[Walter Roberson]

piecewise() is part of the Symbolic Toolbox. It does not accept logical values as its first argument, and does not even accept sym() of a logical value as its first argument (those are converted to numeric values.)

It does accept MuPAD's TRUE or FALSE as its first argument, so if it is important to you to use piecewise() instead of one of the alternatives, then you could use

      FALSETRUE = [sym('FALSE'), sym('TRUE')];

      c = piecewise( FALSETRUE((a2>0)+1), a1^2, 0);

but really you would be better off using just an 'if' or at most logical indexing

      c = 0;
      if a2 > 0
         c = a1.^2;
      end

or

      c = (a2 > 0) .* (a1.^2);              %caution: this is not valid if a1 can be inf

[Rupamathi Jaddivada]

Thanks for the quick reply. My code actually involves two parts: 1. Symbolic manipulation and writing all the expressions into a '.m' file. This '.m' file looks like the one I posted in my question. 2. Running the '.m' file that got created automatically using the command 'eval'

I initially define a variable using the second alternative that you mentioned. All of these expressions get manipulated symbolically and the final expressions get written onto a '.m' file. MATLAB writes these conditional expressions as piecewise functions. When I run this file, I face issues that I mentioned in my previous post.

I also tried using the TRUEFALSE method that you elaborated. It seems like this does not work. MATLAB gives out the following error: Cannot prove '0 < a' literally. To test the statement mathematically, use isAlways.

I also tested using 'isAlways' around the condition, but this always makes the expression false and assigns a value of 0.

Please let me know if there is any efficient way to solve this problem.

Thanks, Rupa

[Walter Roberson]

"Running the '.m' file that got created automatically using the command 'eval'"

Symbolic expressions are not written in MATLAB language, so you should never eval() them. You should use matlabFunction() to create function handles to functions that work numerically, or you should subs() specific values into the symbolic expressions. Or if you have a function that depends upon symbolic functions then make sure you pass symbolic arguments, using sym() to do the conversion if necessary.

[Rupamathi Jaddivada]

Sorry, I think I was not clear enough.

I have a function FunctionABC which does symbolic manipulations and writes the symbolic expressions to FunctionXYZ. It looks something like follows:

function functionA
  syms a1 real 
  ......
  %exp1 and exp2 are complex symbolic expressions found 
  %through the code above
  z = (a1>0)*exp1 + (a1<=0)*exp2

    %code to write text into a file which should then be executed
    fileID = fopen('functionXYZ.m','w');
    %code to write first line of function in the new file
    fwrite(fileID, 'function dx = functionXYZ(x)');
    %code to write rest of the function
    fwrite(fileID,'a1=x(1);\n a2 = x(2);\n');
    fwrite(fileID, ['dx(1) =' ,char(z), '*2;\n']);
    fwrite(fileID, ['dx(2) =' ,char(z), '*4;\n']);
    fwrite(fileID, 'end\n');
    fclose(fileID);

    x0 = randn(2,1);
    eval([strcat('FileNameXYZ,'(x0)']);

end

The 'functionXYZ' that gets created through the code looks something like follows:

    function dx = functionXYZ(x) 
       a1 = x(1);
       a2 = x(2);

       dx(1) = piecewise(a1>0,exp1,exp2) * 2;
       dx(2) = piecewise(a1>0,exp1,exp2) * 4;
    end

In functionABC, the eval command executes functionXYZ with some random x. functionXYZ should be able to compute dx, but it is not able to do that because of the piecewise function that gets created automatically by MATLAB in functionXYZ.

I was wondering if I can make MATLAB process it as an if-alse statement rather than piecewise function.

Thanks, Rupa

[Walter Roberson]

Do not eval().

clear functionXYZ
functionXYZ(x0)

And since you are writing out the piecewise yourself, write it as if or logical indexing. If you really really need piecewise then

functionXYZ( sym(x0))

[Rupamathi Jaddivada]

Hi, Thank you very much. But I actually can't directly use functionXYZ because the name of the function is also dynamic.

I'm not writing out piecewise function myself. I write it either like

if a2>0
  z = exp1;
else
   z = exp2;
end

or like:

    z = (a2>0)*exp1 + (a2<=0)*exp2;

Either of this is written in functionABC.

In both the cases MATLAB automatically creates a piecewise function in functionXYZ as z = piecewise(a2>0,exp1,exp2).

So, when you say logical indexing, do you mean the TRUEFALSE method that you mentioned in your answer previously? To reiterate, that method chooses the index based on the condition being true or false. But the problem here is that the condition is a symbolic expression.

If it is something like TRUEFALSE((a2>0)+1) written in functionABC, MATLAB gives an error message saying that 'a2>0 can't be proved literally. Use isAlways'. The problem is that functionABC has only symbolic expressions. It is only in functionXYZ that the numerical evaluations take place.

Also, this is not the entire problem. functionABC is very complex. After defining z with the complex expression, tt actually extracts symbolic information from many other functions and makes substitutions in the conditional expressions. Finally the manipulated simplified expression is what I want to get copied into functionXYZ for numerical evaluation.

Thanks, Rupa

[Walter Roberson]

MATLAB never turns

if a2>0
  z = exp1;
else
   z = exp2;
end

into piecewise(), even if a2 is symbolic. In the case of symbolic, MATLAB would instead complain that it cannot be proved and to use isAlways .

MATLAB does turn (a2>0)*exp1 + (a2<=0)*exp2 into piecewise() for symbolic a2, and it will refuse to matlabFunction a piecewise() call.

If you have a function that has a piecewise() buried in it and you cannot get rid of the piecewise, then instead of calling the function directly, create

   function_wrapper = @(x) real_function(sym(x))

and you might need to

   function_wrapper = @(x) double( real_function(sym(x)) )

[Rupamathi Jaddivada]

Hello, Thanks again.

I'm sorry. Yes, If-else doesn't work at all in this situation. But when I write the conditional expression as

      (a2>0)*exp1 + (a2<=0)*0

It gives an error saying that inequalities can't be multiplied with zero. Then, I tried just

    (a2>0)*exp1

A piecewise function gets produced by MATLAB. Through the debugger, I see that all the arguments of this piecewise function are logical. I don't know why is this happening and this is why MATLAB produces an error while it is trying to evaluate this expression numerically. It gives an error as follows: Undefined function 'piecewise' for input arguments of type 'logical'.

Thanks, Rupa

[Walter Roberson]

Can you attach your code?

[Rupamathi Jaddivada]

Actually, this is all part of object oriented programming. Each class definition has a method like follows:

   methods(Static)
        function [StateSpace,ControlInputEquations] = PVXiaControlDynamics(this)

            %Inputs____________________________________________
            %entire object
            %_________________________________________________

            %Outputs_________________________________________
            %StateSpace = [ diLddt ; diLqdt];
            %_________________________________________________

            iPVd = this.StateVariables(1);
            iPVq = this.StateVariables(2);
            vDCd = this.StateVariables(3);
            vDCq = this.StateVariables(4);

            vPVd = this.PortInputs(1);
            vPVq = this.PortInputs(2);

            Rf = this.Parameters(1);
            Lf = this.Parameters(2);
            Vs = this.Parameters(2);

            Pref = this.SetPoints(1);
            Vref = this.SetPoints(2);

            psid = this.ControllableInputs(1);
            psiq = this.ControllableInputs(2);

            Gain = this.ControllerGains;

            dphidt = this.RefFrameSpeed;
            wb = 377;

            %averaged switch duty ratio
            Sd = (vPVd+vDCd)/Vs;
      Sq = (vPVq+vDCq)/Vs;   

            %Load Dynamics
            diPVddt = wb*(- Rf*iPVd/Lf + dphidt*iPVq + (Vs*Sd - vPVd)/Lf);
            diPVqdt = wb*(- Rf*iPVq/Lf - dphidt*iPVd + (Vs*Sq - vPVq)/Lf);

            %---- Two layer SM control for PV: Layer 2 -------%
      PPVref = Pref + 1*Rf/(Rf^2+Lf^2)*(vPVd^2 + vPVq^2 - Vref^2);

            %---- Dynamic FBLC in TSS ----% 
            vDCd = Sd*Vs; vDCq = Sq*Vs;
      Pe = vDCd*iPVd + vDCq*iPVq;
      v1 = Gain*(Pe-PPVref-Rf*(iPVd^2+iPVq^2))/2;  % compensate the power loss in the filter
      v2 = Gain*(Pe-PPVref-Rf*(iPVd^2+iPVq^2))/2;

            dvDCddt = psid;
            dvDCqdt = psiq;

            psi1exp = (-vDCd*diPVddt + v1)/(iPVd);
            psi2exp = (-vDCq*diPVqdt + v2)/(iPVq);

            psi1c = (abs(iPVd)>1e-3)*psi1exp;
            psi2c = (abs(iPVq)>1e-3)*psi2exp;

            StateSpace = [ diPVddt ; diPVqdt; dvDCddt; dvDCqdt];
            ControlInputEquations = [psi1c; psi2c];
            StateSpace = simplify(StateSpace);
        end
    end

Then there are few other functions which combine all objects together and do some symbolic manipulations. Finally, there is another function that writes all of these equations into a '.m' file as follows:

function xout = PrintMFileSolveSimulate(PS,FileName,varargin)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%The inputs are 
% PS - Power system object
% FileName - The name of the .m file into which the state space equations should 
% be printed and simulated or solved
%varargin - Serves as intitial guess of solving the system of equations % The outputs are
% xin - numerical equilibrium value if input argument solve = 1
%the result is also stored in a .mat  file
% of the same file name prefixed with 'xout_'
% Usage:
% x0 = randn(15,1);
% xin = PrintMfileSolveSimulate(PS,'SMTLLoad.m',x0);
% First command creates a .m file named 'SMTLLoad_solve.m' first in the same folder
% and find the numerical equilibrium
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Extract the name of the file
FileName1 = strsplit(char(FileName),'.');
FileName2 = FileName1{1};
%Start writing the fucntion with the same name extracted
  FileName = strcat(char(FileName2),'_solve.m');
  fileID = fopen(FileName,'w');
  fprintf(fileID, ['function  x1 = ',char(FileName2),'_solve(x0)\n\n']);

%Text to add parameters path
fprintf(fileID, 'addpath(strcat(char(pwd),''/Parameters/''));\n\n');

%Text to add Load parameters values, set points and controller gains of all the modules
for i=1:numel(PS.Modules)
    a = PS.Modules{1,i}; b = a.StateVariables(1); 
    c = extractAfter(char(b),'_'); Module = c; 
    fprintf(fileID, ['load(''',char(Module),'.mat'',']);
    indx = 0;
    for j = 1:numel(PS.Parameters)
        e = extractAfter(char(PS.Parameters(j)),'_');
        if strcmp(c,e)
            fprintf(fileID, ['''',char(PS.Parameters(j)),''',']);
            indx = indx+1;
            if indx>10
                fprintf(fileID,'...\n');
                indx = 0;
            end
        else
            continue;
        end
    end

....

......
%Print the code to get the inital guess x0
if ~isempty(varargin) && (numel(varargin{1})==TotStates)
    x0 = varargin{1};
else
    try
        load(strcat(strcat('x0_',char(FileName1)),'.mat'));
    catch
        x0 = randn(TotStates,1);
    end
end

%Print the code to solve system of equations - compatible the MATLAB
%2016 only. Might require modifications in options for other versions
    fprintf(fileID, 'fhandle = @SolveDynamics;\n');
    fprintf(fileID,['options = optimoptions(@fsolve,''Display'',''iter'',''SpecifyObjectiveGradient'',false);\n',...
        'options.MaxFunctionEvaluations = 100000000; ',...
        'options.MaxIterations = 1000000;\n',...
        'x = fsolve(fhandle,x0,options);\n\n']);

      fprintf(fileID, 'function dx = SolveDynamics(x)\n');

%Start printing the state space needed to solve the equations
%Split x vector
%State variable
for i=1:numel(PS.StateVariables)
    fprintf(fileID,[char(PS.StateVariables(i)), ' = x(', num2str(i), ');','\n']);

%Print controllable input equations
for i=1:numel(PS.ControllableInputs)
    fprintf(fileID,[char(PS.ControllableInputs(i)),' = ' ,char(PS.ControlInputEquations(i)),';','\n']);
end
....
....
%Print state space equations
for i=1:numel(PS.StateVariableDerivatives)
    fprintf(fileID,[char(PS.StateVariableDerivatives(i)),' = ' ,char(PS.StateSpace(i)),';','\n']);
end

%Return dx vector
fprintf(fileID,'dx = [');
for i=1:numel(PS.StateVariableDerivatives)
    fprintf(fileID,[char(PS.StateVariableDerivatives(i)),'\n']);
end
for i=1:numel(PS.DesiredStateVariableDerivatives)
    fprintf(fileID,[char(PS.DesiredStateVariableDerivatives(i)),'\n']);
end
fprintf(fileID,'];\n');
fprintf(fileID, 'end\n');

%Print the code to save data into a MAT file
fprintf(fileID, ['save(''xout_',char(FileName1(1)),'.mat'');\n']);

%Print last few statement to assign output variable
% and then Run the file that is just printed out through the code above
    fprintf(fileID, 'x1 = x;\n');
    fprintf(fileID, 'end\n');
    fclose(fileID);
    eval(['xout=',strcat(char(FileName2),'_solve'),'(x0);']);
end

The code is too lengthy. But maybe you can just concentrate on the variable 'ControlInputEquations'. In the class defintion that I wrote down initially, I have 'psi1c' and 'psi2c' as conditional expressions. When I look at the variable COntrolInputEquations when this object is instantiated, I see that it ControlInputEquations is an array of two piecewise functions with all its arguments being logical. The first step is to resolve this.

If this is possible to do, the same expression probably gets written in the '.m' file automatically which needs to be solved for numerically.

Please let me know if something is unclear.

Thanks, Rupa