Evaulate Piecewise MATLAB function

Question

Rupamathi Jaddivada 2016-12-31

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/318698-evaulate-piecewise-matlab-function

评论： Rupamathi Jaddivada 2017-1-2

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

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Walter Roberson 2016-12-31

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/318698-evaulate-piecewise-matlab-function#answer_249061

在 MATLAB Online 中打开

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

9 个评论
显示 7更早的评论隐藏 7更早的评论

Rupamathi Jaddivada 2017-1-1

在 MATLAB Online 中打开

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

Rupamathi Jaddivada 2017-1-1

在 MATLAB Online 中打开

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 2017-1-2

Can you attach your code?

Rupamathi Jaddivada 2017-1-2

在 MATLAB Online 中打开

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

请先登录，再进行评论。

Evaulate Piecewise MATLAB function

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

采纳的回答

9 个评论
显示 7更早的评论隐藏 7更早的评论

更多回答（0 个）

另请参阅

类别

标签

Community Treasure Hunt

Evaulate Piecewise MATLAB function

0 个评论 显示 -2更早的评论隐藏 -2更早的评论

采纳的回答

9 个评论 显示 7更早的评论隐藏 7更早的评论

更多回答（0 个）

另请参阅

类别

标签

Community Treasure Hunt

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

9 个评论
显示 7更早的评论隐藏 7更早的评论