Search code examples
matlabplotconstraintsnonlinear-functions

I don't know hot to plot this nonlp's feasible reagion in matlab


I searched and watched how to plot the 3 dimensions of nonlp program's

but I still don't know how to plot these constraints.

x^2+y^2+z^2=1

x>=2*y
2*y>=3*z
x>=3*z

I programmed like this but it doesn't work and I think this is wrong program for upper constraints.

func1 = @(x,y,z) sqrt(1-x.^2-y.^2);
func2 = @(x,y,z) max(x-2*y,0);
func3 = @(x,z) max(x-3*z,0);
ezsurf(func1,[0 1 0 1]);
hold on;
ezsurf(func2,[0 1 0 1]);
hold on;
ezsurf(func3,[0 1 0 1]);
axis([0 1 0 1 0 1]);

Solution

  • A way of doing that is to work with volumetric data.

    The idea is to create a 3D space, with a F value. F will be positive in one side of your main equation and negative in the other.

    [X,Y,Z]=meshgrid(-2:0.05:2,-2:0.05:2,-2:0.05:2);
    
    F=X.^2+Y.^2+Z.^2-1;
    

    Then the conditions are applied to this F. If you want to see the object "cut" then substitute the place where conditions are not met with nan. If you want to see a filled object, then substitute it with a positive number (it may need to be negative in another situation, so check this for other examples.

    % Conditions
    cond1=X>=2*Y;
    cond2=2*Y>=3*Z;
    cond3=X>=3*Z;
    
    % If you want the boundaries to show
    F1=F;
    F1(~cond1)=1;
    F1(~cond2)=1;
    F1(~cond3)=1;
    
    % If you want the boundaries not to show
    F2=F;
    F2(~cond1)=NaN;
    F2(~cond2)=NaN;
    F2(~cond3)=NaN;
    

    Then the idea is to create an isosurface in the zero level (your original function). Here in this code I created 2 plots so you can see the differences between filling and not filling the boundaries. Also some fancy coding in order to get nice figures. However the isosurface and patch parts are relevant to obtain the surface and plot it.

    subplot(121)
    iso1=isosurface(X,Y,Z,F1,0);
    p=patch(iso1);
    isonormals(X,Y,Z,F1,p);
    set(p,'FaceColor','red','EdgeColor','none');
    daspect([1 1 1])
    axis equal
    
    camlight
    lighting gouraud
    
    
    subplot(122)
    iso2=isosurface(X,Y,Z,F2,0);
    
    p=patch(iso2);
    isonormals(X,Y,Z,F2,p);
    set(p,'FaceColor','red','EdgeColor','none');
    axis equal
    
    daspect([1 1 1])
    camlight headlight
    lighting gouraud
    

    Result:

    enter image description here

    Additionally, You can calculate the isosurface of the initial boundary F and plot it with a low alpha value, so you can better understand what piece of the domain you are actually selecting.

    enter image description here