close all
clear all
clc
syms x1 x2 M;              %M为罚因子。
m(1)=1;                    
c=8;                       %c为递增系数。赋初值。
a(1)=0;
b(1)=0;                              
f=x1^2+x2^2+M*(4.6084*x1^5-x2)^2;%外点罚函数
f0(1)=150;

%求偏导、Hessian元素
fx1=diff(f,'x1');
fx2=diff(f,'x2');
fx1x1=diff(fx1,'x1');
fx1x2=diff(fx1,'x2');
fx2x1=diff(fx2,'x1');
fx2x2=diff(fx2,'x2');

%外点法M迭代循环
for k=1:100                                                         
    x1=a(k);x2=b(k);M=m(k);
    %牛顿法求最优值
    for n=1:100                                                      
        f1=subs(fx1); %求解梯度值和Hessian矩阵
        f2=subs(fx2);
        f11=subs(fx1x1);
        f12=subs(fx1x2);
        f21=subs(fx2x1);
        f22=subs(fx2x2);
        if(double(sqrt(f1^2+f2^2))<=1e-6)                         %最优值收敛条件
            a(k+1)=double(x1);b(k+1)=double(x2);f0(k+1)=double(subs(f));
            break;
        else
            X=[x1 x2]'-inv([f11 f12;f21 f22])*[f1 f2]';
            x1=X(1,1);x2=X(2,1);
        end
    end
if(double(sqrt((a(k+1)-a(k))^2+(b(k+1)-b(k))^2))<=1e-6)&&(double(abs((f0(k+1)-f0(k))/f0(k)))<=1e-6)   %罚因子迭代收敛条件
      %输出最优点坐标，罚因子迭代次数，最优值
      a(k+1)   
      b(k+1)
      k
      f0(k+1)
      break;
else
      m(k+1)=c*m(k);
end
end

% 绘制目标函数曲线图
xx1=0:50:600;
xx2=130:20:150;
for i=1:length(xx1)
    for j=1:length(xx2)
        if(4.6084*xx1(i)^5-xx2(j)<=0 )
            Z(i,j)=xx1(i)^2+xx2(j)^2;
        else
            Z(i,j)=0;
        end
    end
end
figure(1);
surf(xx1,xx2,Z);
axis([0 600 130 150 0 600])
title('目标函数曲线图');
xlabel('x1');
ylabel('x2');



运行MATLAB后一直出现这个问题，怎么办呀各位？？

[此贴子已经被作者于2016-4-10 10:24编辑过]

解决了